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The  standard  results  of  linearised  optimal  control  theory  are  explored  and 
examined  to  see  how  they  can  be  applied  to  flight  control  systems.  A  feedback 
control  system  operating  on  the  elevators  of  an  aircraft-like  gliding 
projectile  is  investigated.  The  projectile  is  required  to  closely  follow  a 
predetermined  maximum  range  trajectory  in  the  face  of  initial  disturbances. 

Equations  of  motion  are  set  up  and  linearised.  Approximate  solutions  for  the 
maximum  range  trajectory  are  given  and  approximate  analytical  expressions  for 
the  eigenvalues  of  the  plant  matrix  are  derived.  After  assigning  weighting 
values  to  the  state  and  control  variables  in  the  integral  quadratic 
performance  index,  solutions  to  the  Riccati  matrix  equation  are  computed  and 
used  to  evaluate  optimal  state  feedback  gain  vectors.  The  effect  of  this 
optimal  feedback  on  glider  performance  is  observed  from  computed  trajectory 
simulations.  An  optimal,  feedback  gain  vector  is  selected  subject  to 
limitations  on  angle  of  attack,  elevator  deflection  angle,  and  attitude. 

The  question  of  -ncomplete  state  variable  feedback  is  considered  in  the 
interest  of  simpler  engineering  design.  Using  a  reduced  order  system 
representation,  relationships  between  performance  index  weights  and  closed 
loop  poles  are  established  and  a  sub-optimal  system  based  on  feedback  of  only 
one  state  variable  is  investigated. 
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t.  INTRODUCTION 


1 n  recent  years,  technological  advances  have  opened  'i|>  new  concepts  in  t  lie 
design  ot'  flight  control  systt  is  tor  aircraft,  guided  mi  ss  i  Ion  ,  and  military 
projectiles.  The  products  of  these  advances,  Which  arc  now  manifest  m 
sensors ,  microprocessors,  anil  computing  tcchuii|u<':; ,  are  already  hemp 
exploited  in  modern  advanced  aircraft  ami  puttied  missile  design.  Now,  with 
the  dove  I npment  of  smaller  robust  sensors  together  with  small  and  even  cheaper 
mi  crop rnroKRor  chips,  the  possibility  of  monitoring  and  controlling  the  flight 
of  such  military  projectiles  us  bombs ,  shells,  and  artillery  rockets,  is 
emerging. 

Those  new  developments  have  brought  with  them  now  approaches  in  control  system 
design  which  are  particularly  relevant  to  flight  systems.  In  flight,  when 
speed  and  altitude  vary,  system  parameters  are  often  time-dependent  and  often 
per lornta nee  over  only  a  limited  time  of  flight  is  of  interest.  Modern  control 
system  theory  can  take  account  of  such  time  variations  and  limited  times  of 
flight.  Traditional  classical  control  system  design  theory  is,  in  the  main, 
applicable  to  linear  systems  with  l ime~ invar iant  parameters  and  these  systems 
are  required  to  be  asymptotically  stable  as  time  becomes  large.  The  basic 
mathematical  tool  used  in  t  hi'  traditional  methods  is  the  l.aplacc  transform 
whereas  matrices  arc  basic  to  modern  control  theory.  The  handling  of  mail  ires 
has  been  facilitated  by  the  interactive  computer  aided  design  packages  which 
arc  now  becoming  available  and  which  really  have  turned  t hi'  use  of  modern 
control  theory  hv  engineers  into  a  practical  proposi  t  ion.  *  , •* ,  Such  matrix 
methods  can  readily  cope  with  the  mul t i - input -out  put  control  systems 
tieipieully  needed  in  flight,  as  for  instance  on  an  aircraft  with  rudder, 
elevator,  and  aileron  controls.  For  systems  which  arc  sensibly  time- 
invariant,  matrix  methods  can  be  readily  reduced  to  the  well  known  classnal 
representation  thus  rendering  such  systems  amenable  to  treatment  by  the  well- 
developed  classical  design  techniques. 


Another  attraction  of  the  modern  theory  in  applications  to  flight  control 
system  design  is  that  it  has  a  natural  connection  with  optimisation  theory, 


and  this  in  turn  enables  characteristics  of  the  desired  flight  path  to  be  used 
directly  as  criteria  while  developing  a  feedback  control  design.  For  example, 
in  the  case  of  a  feedback  control  system  for  a  projectile  which  is  required  to 
follow  a  desired  trajectory,  the  modern  method  could  be  used  to  find  a  set  of 
optimal  feedback  gains  which  would  serve  to  minimise  the  time  integral  of 
weighted  squares  of  deviations  from  the  required  trajectory  values.  Another 
use  would  be  to  determine  optimal  feedback  gains  for  a  terminal  guidance 
system  to  minimise  miss  distance  subject  to  limitations  on  the  amount  of 
available  control.  A  virtue  of  the  linear  control  system,  which  results  from 
this  method  of  minimising  the  integral  of  weighted  squares  of  errors  of 
deviations  in  system  variables,  is  that  by  following  simple  rules  governing 
selection  of  weights  the  system  is  inherently  stable  and  also,  the  effects  of 
random  variations  or  noise  are  suppressed.  In  fact,  the  art  of  linear  system 
design  by  this  method  often  lies  in  the  choice  of  suitable  weights.  In 
general,  the  optimal  feedback  gains  will  he  time-varying  and  an  on-board 
microprocessor  will  be  needed  to  generate  the  gain  values  together  with 
desired  values  of  system  variables  to  be  compared  with  actual  values  measured 
by  sensors. 

This  thesis  aims  at  exploiting  modern  control  theory,  to  investigate  how  it 
ran  be  applied  to  the  mid-course  flight  of  a  gliding  projectile  and  to  see  how 
it  reveals  the  physical  flight  behaviour.  The  sort  of  projectile  considered 
is  a  standard  bomb  configuration  to  which  an  aircraft-like  wing  has  been  added 
together  with  a  pair  of  elevator  controls  at  the  rear  tail  fins  as  shown  in 
figure  1.  The  wing  provides  the  bomb  with  a  gliding  capability,  significantly 
increasing  its  range  so  that  the  position  of  launch  may  remain  at  a  safe 
stand-off  distance  from  a  heavily  defended  target.  A  complete  trajectory  from 
launch  to  target  would  consist  of  three  phases,  as  shown  in  figure  2,  launch, 
mid-course,  and  terminal.  Launch  may  occur  by  release  from  an  aircraft  or  by 
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rocket  boost  from  a  ground  launcher  and  the  launch  phase  of  the  trajectory 
must  terminate  in  favourable  starting  conditions  for  the  mid-course  gliding 
flight.  Mid-course  flight  must  take  the  projectile  to  within  acquisition 
range  of  the  target,  such  that  the  terminal  phase  seeker  can  see  the  target. 
However,  ttie  greater  the  mid-course  range  the  greater  the  terminal  errors  due 
to  initial  disturbances  will  be,  and  hence  the  need  for  some  corrective 
control,  and  for  some  assessment  of  allowable  tolerance  in  mid-course  initial 
conditions . 

The  main  influence  of  elevator  controls  is  on  pitching  motion  which,  for  an 
aircraft  at  zero  bank  or  roll  angle,  means  that  the  main  influence  is  on  those 
system  variables  which  affect  motion  in  a  vertical  plane.  Such  motion  is 
termed  longitudinal  motion  as  distinct  from  the  lateral  motion  which  governs 
roiling  and  yawing  behaviour.  For  conventional  aircraft  it  is  well- 
established  that  the  longitudinal  and  lateral  motions  can  be  decoupled  as  long 
as  roll  angle  variations  do  not  exceed  first  order  magnitudes.4,5  On  the 
ideal  mid-course  trajectory  where  there  is  no  call  for  side  manoeuvre  there  is 
no  demand  for  large  roll  angles  and  hence  for  the  purpose  of  preliminary 
control  system  design,  the  lateral  and  longitudinal  motions  can  be  treated 
separately.  For  aircraft,  the  longitudinal  motion  typically  possesses  two 
characteristic  modes,  a  well-damped  short  period  mode  and  a  lightly  damped 
long  period  mode.  Oscillitary  motion  of  the  short  period  mode  is  prominent  in 
the  angle  of  attack  and  angular  pitching  velocity  behaviour  whilst  the  long 
period  oscillations  dominate  the  velocity  and  attitude  behaviour.  The  damping 
agency  for  the  long  period  mode,  or  phugoid  as  it  is  called,  is  the 
aerodynamic  drag.  In  order  to  maximise  range  it  is  necessary  to  minimise  drag 
and  hence,  damping  is  of  necessity,  small.  This  means  that  initial  errors  in 
velocity  and  attitude  would  tend  to  persist  throughout  mid-course  flight 
without  some  form  of  control.  Using  figure  3,  a  simple  physical  explanation 
of  the  phugoid  oscillation  can  be  given.  For  the  undisturbed  or  equilibrium 
flight  path,  the  vertical  component  of  aerodynamic  lift  is  a  little  less  than 
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the  weight,  allowing  the  glider  to  sink  at  a  rate  at  which  the  potential 
energy  loss  just  balances  the  energy  dissi paterl  by  drag.  Suppose  that  the 
glider  flies  at  a  constant  angle  of  attack,  then  a  disturbance  which  causes 
the  velocity  to  increase  will  cause  the  lift  to  increase  and  the  glider  will 
rise,  increasing  its  potential  energy  at  the  expense  of  its  kinetic  energy. 
Velocity  then  decreases,  causing  lift  to  decrease  until  the  glider  eventually 
sinks  and  the  process  is  repeated  in  reverse  setting  up  a  phugoid  oscillation 
about  the  equilibrium  flight  path.  Aerodynamic  drag  dampens  this  oscillation. 
To  counter  this  oscillation,  the  angle  of  attack  would  need  to  be  continuously 
controlled  so  that  the  equilibrium  balance  between  lift  and  weight  of  the 
undisturbed  flight  path  is  maintained.  This  is  what  the  elevator  controls 
seek  to  do  by  their  pitching  effect  on  angle  of  attack. 

In  a  powered  aircraft  the  phugoid  would  be  controlled  by  regulating  velocity 
through  engine  throttle  control,  relying  on  elevators  for  attitude  control. 
In  a  piloted  glider  the  pilot  is  usually  able  to  counter  phugoid  motion 
because  of  its  relatively  long  period.  The  unpiloted  glider  relying  on 
automatic  control  of  phugoid  motion  through  elevators  has  not  attracted  much 
attention  in  the  literature  although  it  is  interesting  to  note  that  automatic 
control  systems  which  fed  back  attitude  and  velocity  to  elevator  controls  were 
devised  for  both  British  and  German  bomber  aircraft  early  in  the  second  world 
war.6  When  a  gliding  projectile  is  to  be  released  from  an  aircraft  it  is 
likely  that  the  error  in  initial  attitude  reference  will  exceed  the  magnitude 
of  attitude  measurements  needed  for  feedback  during  gliding  flight,  so  it  is 
worthwhile  briefly  considering  the  consequences  of  sensing  velocity  alone. 
One  way  of  counteracting  the  phugoid  is  for  the  elevators  to  cause  the  glider 
to  climb  for  velocities  in  excess  of  equilibrium  and  attempt  to  maintain  the 
increased  altitude,  when  velocity  falls  to  the  equilibrium  value.  Although 
this  diminishes  the  phugoid  results  in  a  displaced  equilibrium  flight,  path 
and  loses  positional  act  :.-cy.  To  maintain  positional  accuracy  it  is 
necessary  to  try  to  inaiii'-.iu  oscillations  about  the  equilibrium  flight  path. 
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This  could  lx*  achieved  by  allowing  angle  of  attack  to  dec. ease  for  velocities 
in  excess  of  equilibrium  values,  thus  diminishing  the  excess  of  lift,  over 
weight  and  hence  phugoid  amplitude,  but.  also  diminishing  drag  and  hence 
damping.  This  tends  to  maintain  equilibrium  flight  path  position,  but  unless 
the  control  system  is  carefully  designed,  too  large  a  decrease  in  angle  of 
attack  in  response  to  increased  velocity  could  lead  to  an  unrecoverable 
downward  plunge.  In  this  respect,  phugoid  motion  is  not  unlike  that  of  a 
simple  pendulum  which  executes  small  oscillations  about  its  equilibrium 
position  for  small  disturbances,  but  for  large  disturbances  a  projectile  can 
loop  the  loop  whereas  the  pendulum  encircles  the  pivot. 


! 
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If,  on  the  other  hand,  the  angle  of  attack  is  caused  to  increase  for 
velocities  greater  than  the  equilibrium  value  then  the  opposite  effects  occur 
and  phugoid  oscillations  are  enhanced.  The  amplitude  of  oscillations  is 
initially  increased  but  the  increased  angle  >  attack  increases  the  drag  and 
hence  the  damping  has  a  more  marked  effect  in  reducing  excess  velocity. 

The  consequences  of  the  foregoing  discussion  arc  that  we  can  expect  to 
regulate  velocity  and  attitude  at  the  expense  of  position,  velocity  and 
position  at  the  expense  of  attitude,  or  attitude  and  position  at  the  expense 
of  velocity.  These  observations  will  be  useful  later  on  when  weights  are 
selected  for  constructing  the  performance  indexes  used  in  optimal  feedback 
design.  For  the  projectile  considered  in  this  thesis,  the  objects  of  the  mid¬ 
course  control  then  are  to  maximise  range  by  preserving  gliding  flight  and,  to 
within  constraints  on  control  deflection  and  angle  of  attack,  to  reduce 
terminal  effects  arising  from  initial  errors  and  flight  disturbances. 

Elements  of  linearised  optimal  control  theory  are  introduced  in  sect. on  5  and 
application  of  this  theory  to  the  design  of  a  feedback  control  system  for  the 
glider  is  carried  out  in  section  6.  A  consequence  of  linearised  optimal 
control  theory  is  that  all  state  variables  should  be  fed  back  to  the 
controller  thus  it  is  of  interest  to  investigate  the  performance  of  simpler 
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sub  optima]  systems  wi th  partial  feedback.  A  sub  optimal  feedback  control 
system  based  on  a  reduced  order  mathematical  model  which  considers  only 
phugoid  motion  in  isolation  is  studied  in  detail  in  section  7.  Also 
considered  in  section  7  is  a  sub  optimal  system  in  which  only  velocity  is  fed 
back  to  the  controller. 

In  the  next  section  the  equations  of  motion  for  a  gliding  projectile  are 
derived  and  include  the  effects  of  wind  disturbances.  In  section  3,  the 
equations  of  motion  are  used  to  extract  conditions  for  optimal  flight  paths 
and  in  section  4  the  equations  for  first,  order  perturbations  about  optimal 
conditions  are  derived  and  form  the  basic  equations  for  applying  the  linear 
optimal  control  theory  in  the  section  following  on. 

2.  EQUATIONS  OF  MOTION 

1  he  full  equations  of  Motion  tor  an  aircraft  can  be  obtained  from  sncli  texts 
as  references  4  or  5 .  It  is  found  that  when  large  sideways  manoeuvres  are  not 
needed  then  roll  angles  can  remain  small,  and  horizontal  and  vertical 
components  of  the  motion  become  uncoupled.  Thus  for  the  purpose  of 

preliminary  design  of  a  feedback  control  system  operating  on  the  projectile's 
elevators,  only  motion  in  a  vertical  plane  will  be  considered.  Equations  of 
motion  for  flight  in  a  vertical  plane,  subject  to  wind  disturbances,  will  now 
be  derived  in  a  form  convenient  for  treating  gliding  flight  and  for  subsequent 
linearisation  for  the  purpose  of  control  system  design.  Nomenclature  together 
with  force,  moment,  and  axes  systems  for  flight  in  a  vertical  plane  are  shown 
in  figure  4.  In  the  figure 

V  is  the  velocity  vector  of  the  projectile's  centre  of  gravity, 

Vy  is  the  velocity  vector  of  the  atmospheric  wind, 
and  =  -V+Vy,  the  velocity  vector  of  the  air  relative  to  the 

projectile  and  thus  the  velocity  on  which  the  magnitude  and 
direction  of  the  aerodynamic  forces  depend. 


The:  aerodynamic  forces  a  re :  - 


1.  the  lift,  force  acting  normal  to  \f  not  including  tlio  elevators 
component. , 

1..,  the  lift  force  component,  of  the  elevators, 
and  l)  the  drag  force  acting  in  the  direction  of  V^. 

N  is  the  pitching  moment  about  the  centre  of  gravity, 
and  m  is  the  mass  of  the  projectile. 


The  angles  in  figure  4  are:- 


« 


Of 


R 


•1 

and 


the  angle  of  attack,  which  is  the  angle  between  V  and  the 
projectile's  longitudinal  axis,  positive  for  nose  up, 
the  re'ativo  angle  of  attack,  which  is  the  angle  between  and  the 
projectile’s  longitudinal  axis, 

the  elevator  control  deflection  angle,  positive  downwards  as  shown 
8  the  angle  between  the  projectile's  longitudinal  axis  and  the 
horizontal  lying  in  the  same  vertical  plane. 


There  are  two  sets  of  three  dimensional  axes  systems  shown  in  figure  4  and 
these  are  defined  by  the  mutually  perpendicular  unit  vector  triads  £  , 
and  C ,  j ,  n . 

The  system  is  fixed  relative  to  earth  with  ^  and  j*  contained  in  a 

horizontal  plane  such  that  is  directed  along  the  range  X  and  is 
vertically  downwards.  The  t,  j,  n  system  with  origin  at  the  centre  of 
gravity,  rotates  with  angular  velocity  j (0-rt)  and  t  is  along  the  tangent  to 
the  flight,  path  and  n  is  along  the  inward  drawn  normal. 


The  linear  acceleration  of  the  centre  of  gravity  is  found  by  differentiating 
the  velocity  vector  V  =  Vt . 


*  •  • 

Thus  V  =  Vt  +  Vt. 

where  t  ~  j(ft-of)  AC  =  -h(^-ot), 


giving  t!u>  result  V  -  V(  -  V(t>*«)  a. 

From  figure  4,  the  force  acting  on  the  projectile  is 

-  {!)  cos  -  (1,+ly)  sin  ♦  mg  sin  ((-)-«) }  f 

-  {(L+h^,)  cos  (cy^-w)  +  D  siii(a^-ot)  -  mg  cos  (O-a)jii 

» 

Equating  this  force  to  in^-  gives  the  two  scalar  equations 

* 

mV  =  -Dcos(tf„-«)+(L+L  )sin(aB-0()-mg  sin (0-tv)  (1) 

K  C  K 

4  « 

ami  mV(0-a)  =  (I.+L  )cos(aD-«)+Dsin(«  -«)-ine  cos(,0-a)  (2) 

C  K  K 

In  addition  there  is  the  moment  equation 

B0  =  M  (3) 

where  B  is  the  pitching  moment  of  inertia  of  the  projectile  about  a  transverse 
axis  through  its  centre  of  gravity. 

Additionally,  the  range  X  and  altitude  h  satisfy 


X  =  V  cos 

(0-«) 

(4) 

h  “  V  sin 

(0-«) 

(5) 

It  is  now  necessary  to  define  the  aerodynamic  forces  and  moments  appearing  in 
equations  (1),  (2)  and  (3).  The  drag  is  expressed  in  the  forn/^ 

D  =  (pV*S/2) 

x{Cl)o+0fKS(3lV30,)2/(”KASw)} 
is  air  density 

is  the  reference  area  on  which  aerodynamic  coefficients  are 


where  p 
S 

based , 

C 


no 


is  the  drag  coefficient  at  zero  angle  of  attack, 


OC j. / Oof  is  the  wing  lift  coefficient  derivative  evaluated  at  tt~Q 

and  S/(;ikAJ?  )  is  a  constant, 

w 

The  «R  term  in  the  expression  for  drag  1),  represents  the  component  of  drag 
induced  by  vort.ici  t.y  shed  from  a  finite  span  wing  and  would  disappear  only  for 
a  wing  of  infinite  span.  The  lift  on  the  projectile,  apart  from  that  arising 
from  deflected  elevators,  is  given  by 

1,  =  (pVj*S/2)  CL(«R). 

and  the  lift  on  the  deflected  elevator  control  by 

t,c  =  (ijv*s/2)  (ocL/&n)n, 

where  CT («u)  is  the  lift  coefficient  for  the  projectile 

and  c)CT /an  is  the  lift  coefficient  derivative  for  the  elevator  control 
evaluated  for  t]-0. 

The  drag  coefficient  C  and  the  lift  coefficient  derivatives  f>C,  /Oij  and 
DCj/9«  are  Mach  number  dependent,  but  for  the  Mach  number  range  considered 
here,  which  is  less  than  O.T),  these  aerodynamic  coefficients  can  reasonably  be 
taken  as  constants.  Thus  the  projectile  lift  will  be  expressed  as 

(pvJs/2)  Idc,/0«)«R 

where  (OCj/Ort)  is  evaluated  at  Of~0 . 

Similarly,  the  pitching  moment  is  expressed  as 

M  =  (pVgSb/2)  {«R(0t:rn/;v.)!(qh/2VR)(3cm/i)«i) 

+(a  b/2V  )(00  /<)n)+n(c)C  /i)n)} 

K  h  m  m 

where  b  is  a  reference  length 

C  is  the  aerodynamic  pitching  moment  coefficient 
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and  q  -  0 

and  the  partial  derivatives  of  C  are  evaluated  for  zero  values  of  their 
r  m 

arguments.  The  moment  components  in  0C^/D«  and  3C  /3ij  represent  the  static 
pitching  moment  due  to  angle  of  attack  and  elevator  deflection  respectively, 
while  the  components  ill  9C  /<)q  and  represent  damping  effects. 

The  air  density  p  appearing  iti  the  aerodynamic  torees  and  moments  is  dependent 
upon  altitude  h  and  a  frequently  used  analytical  relationship  is  given  hy 

-h/H  . 

P  =  Pr<>1-  e  ,  where 

p  ..is  a  constant  reference  density 

re  t 

and  li  is  a  constant  scale  height. 

Using  this  relationship  for  p  and  expressing  the  relative  velocity  VR  and 
angle  of  attack  oR  in  terms  of  wind  velocity  Vy  and  wind  direction  with  the 
horizontal  p^,  as  shown  in  figure  /«,  equations  (1)  to  (T>)  together  with 

0  -  q 

give  the  following  set  of  equations  for  the  motion  of  the  projectile. 

V  =  -(D/m)  {V-Vw  cos(M+0-«)}  /  VR 

-(I./m)(V  /VD )  sin(|*  +0-n)  -g  sin(O-cr) 

W  K  W 

*(p,  fS/2m)(ac  /9q)r|V  V  sin(p  +0-or)exp(-h/H)  (6) 

t  C I  I*  f\ 

«  =  (D/mV)(V-w/VR)  sin(pw*0-rt) 

-(I. /mV)  {V-Vw  cos(»iw»0-u)}  /  VR 
+(g/V)  cos (0-0 )  *  q 

-(prpfs/2n)OcL/an)  nvR{v-vw  cos(p  +o-«)} 


(7) 


n  - 


<*>  =  4 

q  -  (|»  rKb/2B)V£{«tAOC  /0«)»(*|I»/2V»  /<)q) 

i  t  i  ^  **  hi  “in 

» (a-h/m'rtHoc  /c>f»)*n(c)c  /oqi  !e'h/" 

*  *  in'  m 

X  =  V  cos(O-u) 

# 

h  •*  V  sin((-V-tV) 


(B) 

(')) 

(  10) 

(11) 


liquations  (<>)  to  (11)  thus  represent  .i  sot  of  six  equations  for  the  state 
variables  V.  I'i,  0.  q,  X,  h  and  control  variable  I)  where 


«R  =  rt{V;’-VVweost|iw.t-)~u)(  /V" 

►  tM  tq){V-’-W  cos  ( |  i  »0-o)i  /Vj; 
www  w  K 

+  {(VV  -V  V)sin(|J  +0-<v)(  /\’l  (10) 

W  W  W  t\ 


The  values  ot  the  aerodynamic  coefficients  appearing  in  equations  ((>)  to  (13) 
can  be  found  from  data  handbooks  such  as  references  7  and  8. 


3.  OPTIMAL  GLIDK  PATH 


The  desired  glide  path  is  one  which  maximises  the  range  X  for  a  given  decrease 
in  altitude  for  flight  in  still  air.  Descent'  rate  is  governed  by  the  elevator 
control  setting  so  that  the  mathematical  problem  is  to  find  the  time  variation 
of  control  angle  r)(t)  which  maximises 


t. 

,i 


subject  to  the  constraints  imposed  by  equations  (6)  to  (16)  with  0.  When 
Vw= 0,  equations  (6)  to  (16)  become, 


V  =  -(p  ,S/2m)V2Cn(a)e"h/H  -  gsin(e-a)  (17) 

re  r  1J 

a  =  -(prefS/2m)V{<*OCL/9oi) 

+  n(<5CL/8n)}e_h/H  +  o i  +  (g/V)cos  (0-ot)  (18) 

0  =  q  (19) 

q  =  (p  ,Sb/2B)V2  e~h/H{a(9C  /8a)+(qb/2V) (dZ  /8q) 
ret  m  m 

+(&b/2v)(ac  /aa )+oOc  /on)}  (20) 

m  m 


where 

Cu(a)  =  CDq  +  a2S(3Ct/3a)2/At.ASw  (2!) 

However,  a  good  engineering  approximate  solution  is  available  for  this  problem 
and  is  based  on  physical  reasoning.  It  is  argued  that  maximum  range  flight 
paths  will  be  flat  and  hence  changes  in  0  and  a  will  be  small.  Hence,  the 


JVCW  3 3VS3  iV^JJTWtVWSf^fi  V* 
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approximation  G  =  q  =  a  ~  0  is  attempted  in  equations  (17)  to  (21),  giving 


V  -  -(prejrS/2in)V2C[)(«)e  h^H  -g  sin(G-a) 


(22) 


0  =  -  (P,:f,(S/2ni)V{a(DCj  /3a) 


-h/H 


+  n(»CL/3n)}e  +(g/V)  cos (0-a) 


(23) 


0  =  «(3C  /3a)  +  n(3C  /3r|) 

<H  1)1 


(24) 


The  physical  interpretation  of  these  equations  is,  first,  from  (23)  that  the 
lift  is  approximately  equal  to  the  weight  of  the  projectile  and  that  changes 
in  velocity  are  due  to  changes  in  altitude.  For  small  changes  in  velocity,  or 
that  is,  kinetic  energy,  equation  (22)  implies  that  the  rate  of  loss  in 
potential  energy  is  approximately  balanced  by  the  energy  dissipated  in  drag. 
Equation  (24)  is  the  trim  equation  giving  the  relative  values  of  a  and  r|  which 
reduce  the  pitching  moment  to  zero.  As  shown  in  reference  (9),  equations  (22) 
and  (23)  can  be  combined  to  give 


d(V2+gh)/dt  -  -gXCJ>(a)/ta(3CL/3a)+nOCu/;n)} 


where,  from  equation  (24),  f|  -  -a(3Cm/3a)/(0Cn|/3ij) , 


so  that  integration  of  equation  (25)  with  a  assumed  to  be  sensibly  constant 
yields 


X-X.  ■-  { (V2-V2)/2g  +(h.-h)H«/c1)(«)} 
x{  (c)CJ/3a)+(n/a)(aC] /On)  1 


(26) 


where  subscript  i  denotes  initial  value.  It  is  possible  to  obtain  further 
close  approximations  for  trajectory  parameters  wher,  a  is  assumed  to  be 


invariant.  For  example,  equations  (23)  and  (24)  give 


V2  =  g  cos  (<•)-«) 

_  (27) 

<Pr„t-S/ 2m )  { ( 0<V ;>« )  +  ( n/« )  (  <7 Cl  /!) i] )  jo 

Differentiating  equation  (27)  with  respect  to  t i me  and  ignoring  changes  in  ((-)- 
a)  gives 


V  -  ( V‘:/2H )  sin (0-«)  (28) 

which  when  substituted  into  equation  (22)  gives  for  the  glide  path  angle*. 


tan(O-o)  -  -g  CD(«) 

{ g  t V2  /  2H )  { ( 0CL  /  Xi)  +  ( I) /  « )  ( <>Ct/  Oi  | ) }  (29) 

The  time  of  flight  can  be  found  by  putting 


rlh/h 


dh/;V  siii(O-o) 


h 

i 


h  . 
x 


Equations  (27)  and  (29)  enable  (V  sin(O-o)}  ^  to  be  expressed  in  simple 
exponential  functions  of  h  which  can  be  readily  integrated  to  give  the  result 


t  -  {(OCL/ta)  +  (n/n)(dCL/Oi])  } 

x  ojv'V^VJ/g  -2H/V.+2H/V}  /  CD(o)  cos(O-e)  (30) 

Equations  (26)  to  (30)  provide  close  approximations  to  glide  path  trajectories 
which  are  flown  with  angle  of  attack  a  sensibly  constant  such  that  the  glide 
angle  is  not  too  steep  in  the  sense  that  the  approximation  cos(()-nf)  s  1  is 
acceptabi e . 


It  can  be  seen  from  equation  (26)  that  if  changes  in  V‘:  are  ignored  a  maximum 


value  for  the  range  (X  -  Xp  occurs  when  a  constant  angle  of  attack  is  chosen 
to  maximise  «/Cp(o).  ignoring  changes  in  V2  results  in  an  optimum  angle  of 
attack  which  is  independent,  of  the  altitude  range.  From  equation  (2!)  it  is 
found  that  the  if  is  t lie  value  of  a  which  maximises  a/C^(ot),  then  the 
optimal  angle  of  attack  is  given  by 


(S/rt£ASw)(OCL/aa)‘ 


The  foregoing  engineering  approach  shows  that  a  close  approximation  to  the 
glide  path  with  maximum  range  occurs  when  the  angle  of  attack  is  held  at  a 
constant  value  a and  hence,  from  equation  (24),  the  corresponding  constant 
elevator  setting  is 


n  =  -a  (3C  /3a)  /  OC  /3n) 

o  o  m  n 


The  optimum  angle  of  attack  is  a  function  of  the  aerodynamic  properties  of  the 
projectile  and  equations  (27)  and  (29)  show  that  the  initial  optimal  velocity 
and  glide  path  angle  are  dependent  on  aQ  and  initial  altitude  h^ .  Thus  for  a 
selected  value  of  h^,  if  the  launch  phase  cannot  consistently  provide  an 
initial  mid-course  velocity  as  high  as  the  required  optimal  value,  then  ct 
would  have  to  be  increased  from  the  optimum  to  satisfy  equation  (27),  and  in 
consequence  the  range  would  be  reduced.  On  the  other  hand  if  the  launch  phase 
could  consistently  provide  an  initial  mid-course  velocity  in  excess  of  the 
optimum,  decreasing  a  below  the  optimum  would  again  lead  to  reduced  range 
provided  the  assumption  that  changes  in  V2  were  small  still  held.  In  this 
case  som2  other  strategy  would  be  employed  to  retain  potential  energy  until 
sufficient  kinetic  energy  had  been  dissipated  to  slow  the  projectile  down  to 
the  optimal  velocity  and  then  initiate  the  glide.  It  will  be  assumed  that  the 
launch  phase  can  be  designed  such  that  the  expected  conditions  at  the 


commencement  of  mid-course  flight  arc  the  optimal  glide  path  values.  The 


L#>  - 


function  of  the  feedback  control  system  actuating  the  elevators  is  then  to 
minimise  errors  in  position,  at  the  end  of  mid-course  flight,  due  to  initial 
departures  from  optimal  glide  path  trajectory  parameters. 

4.  PERTURBED  EQUATIONS  Of  MOTION 

Departures  from  the  optimal  glide  path  parameters  are  now  treated  as  first 
order  perturbations  about  optimal  values.  A  set  of  linear  equations  in  these 
perturbed  state  variables  can  now  be  found  from  equations  (6)  to  (11).  The 
procedure  is  illustrated  by  taking  equation  (6)  as  an  example.  This  equation 
is  of  the  form 

Fy(V,  a ,  0,  h,  V,  >1,  Vw)  -  0, 

Therefore  6F v  =  0  , 

i.e.  (UFy/DV^V,  +(OFV/Uot)o«1 
+(3FV/U0)o  0J+(9Fv/3h)oh1 

+  (3Fy/yv)0vJ  +(0Fv/an)onJ  +OFv/ov  )vw  -  o  (33) 

where  subscript.  0  means  that  partial  derivatives  are  evaluated  at  optimal 
values.  Carrying  out  this  procedure  for  each  equation  from  (6)  to  (11)  yields 
a  system  of  linear  first  order  equations  in 

*1;  D .  > 

which  are  the  first  order  perturbations  in 

V,  a,  G,  G,  h,  n,  respectively. 

The  resulting  perturbation  equations  are  now  given 


V,  -----  -(p  ,.S/2m)  2V  C.,(n  )V. 

1  ' rel  o  D'  o  l 


+  Of  gcos(0  “Of  ) 
i  o  o 


ai  (prefS/2B)  V2o(aCD/3a)c 


0,  g  cos (0^-ti^) 


+  (prefS/2ro>  (VVH)  S«V  hl 

+  V°S  Pw  (prefS/2ra)  2V«CD(<*o)cos^o-<Yo) 
+Vwco«Mw(prefS/2.)  Vosin(0o-c.o) 

x{OC0/3a)o  -cL(ao)  -n0Oct/an)} 
-Vw«inMw(Prefs/2m)  2VoCD(«o)sin(Oo-oo) 

+VwSinMw(Pre£s/2m)  V()COS  (6^) 

x{(dcD/da)o  -cL (otQ)  -  n0Oct/an)} 


-(gV./V2  )  cos (0  -a  ) 

I  O  O  o 

■Vl(prefS/2m)  {CL(ofo)+flo(aCL/3n)} 

-(a,/Vo)  (PrefS/2m)  V2q  OC^/Ba) 

♦<VV«  sin(©  -or  ) 

-(^g/Vo)  sin(eo-ao)  +  q 

+  (PrefS/2m)  (Vohj/H)  (CL(oo)  -*-n0(3CL/Bn)} 
'(prefS/2ra)  Vo(9CL/3,l)  ni 


+2Vwcospw(prefS/2m)  cos(0o-ao) 

x{cL(ao)  +nQOcL/9n)} 


+V.  c°s|.i  (p^e<.S/2ni)  sin(0_-a  ) 


x{0CL/aa)o  +CD/(tfo)} 


"2VinMw(prefS/2m)  sin(V°'o) 

x  { Cr  (of  )  +n  (3C  /Dr|)} 

L  O  O  Li 


+VinPw(prefS/2m)  COs(eo"°'o) 


[{0CJ/8a)o  +CD(ffo)} 


-vl<PcefSI,2/4DVo)(8Cm/8i,)o 

x(prcls/2"')''02fcL<«0)*n0(acl/3n)) 

-v,  oc „/3i)0 

x(p  _Sb2/4BV  )  g  cost©  -a  ) 
rret  o  oo 

+«.  (p  Sh/2B)V2  (3C  /9a) 

1  vr  ref  o  m  o 

-UPretS/2m)v2o  0CL/9a)o  -g  sin(Oo-ao)} 
xal  (PreiSb2/AB)(3Cm/3,i)o 

-0(p  Sb2/4B)(3C  /9a) 

*  rrei  m  o 

xg  sin(0  -a  ) 
o  o 


+<l(prefSb2/4B)Vo 

x{Ocm/3q)0  HBCJda) o] 


-hl(prefSb2/4B)  ^8Vai)o 
x{prefS/2m)  (Vq/H)  (Cj  («o)+noOCL/9n)) ) 


+n,(p  fSb/2B)V2  (DC  /9n) 

1  ref  o  m 
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+Vw  cosMw(pr(|fSb2/4H) 

xOCw/a&)ocvo/vo)sin(eo-ao) 

-Vw  cosMw(Prel.Sb/2B)V0  (acm/a«)  sin(0  -a  ) 

+V  si up  (p  Sb2/4B) 
w  fw  Kref 

x(OC  /da)  (V  /V  ) cos (0  -a  ) 
m  o  o'  o  o  o 

-vw  si“Hw(PrefSb/2B) 

xV  (9C  /da)  cos (0  -a  ) 

O  111  o  o  o 

-(PrefSb2/4B)  (DCm/9a)o 

xsin(0  -a  )  d/dt(V  cosu  ) 
o  o'  w  rw 

-(PcefSt.V4B)  (3Cm/ai)o 

xcos(eo-«o)  d/dt(Vwsinpw)  (37) 


h,  =  V . sin(0  -or  )  -a  V  cos(0  -a  ) 

1  1  o  o  1  o  o  o 

+6.V  cos(0  -a  )  (38) 

«  o  o  o  v  J 

Equations  (34)  to  (38)  are  in  a  form  amenable  to  application  of  linear 
optimal  control  theory. 


5.  LINEAR  OPTIMAL  CONTROL  THEORY 
Equations  (34)  to  (38)  can  be  written  in  the  standard  form 


x  =  Ax  +  bf|j  +  Dw 


where 


X  is  a  5-  component  state  variable  vector 
A  is  a  5x5  system  matrix 
b  is  a  5-  component  control  vector 

is  increment  in  elevator  control  angle 


(39) 
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1)  is  a  5x5  disturbance  matrix 
and  w  is  a  5-  component  wind  vector 

In  standard  notation,  x  is  the  state  variable  column  vector  having  components 
Xj ,  where  corresponding  to  equations  (34)  to  (38),  j=l,2,3,4,5,  such  that 

x,  =  V,  ,  xa  =  Of,  ,  xA  ^  0r  x„  s  q,  xf  H  h, 

Optimal  control  theory  is  thoroughly  covered  in  many  texts.10,11,16.  In 
reference  10  optimal  feedback  control  is  introduced  through  the  Hamilton- 
Jacobi  approach  which  is  quite  general  and  not  restricted  to  linear  systems. 
Theory  of  linear  optimal  systems  is  treated  in  great  detail  in  reference  11 
while  reference  16  deals  with  optimisation  generally.  Because  this  thesis  is 
restricted  to  linear  optimal  control  of  deterministic  systems  it  is  worthwhile 
demonstrating  how  readily  solutions  for  such  linear  systems  follow  from  the 
standard  results  of  the  calculus  of  variations.  This  leads  directly  to  the 
optimal  system  description  with  the  additional  adjoint  variables  and  to  the 
alternative  Riccati  matrix  equation.  In  the  absence  of  random  disturbances 
the  deterministic  problem  is  treated  as  follows.  First  an  integral  quadratic 
performance  index  of  the  form 

tf 

J(tf,t.)  -  j  (xTQx+n^I>n1)<lt+xT(tf)Sx(tf)  (40) 

t . 
i 

is  set  up  where  Q,  P,  S  are  symmetric  weighting  matrices  which  are  arbitrarily 
selected  hut  subject  to  the  conditions  that  P  is  position  definite  and  Q  and  S 
are  non-negative  definite.  The  object  is  to  find  the  functional  form  of  the 
control  angle  increment  r|j  which  minimises  J.  This  is,  in  effect,  to 
determine  a  control  setting  which  minimises  the  variance  of  perturbations  in 
the  state  variables,  subject  to  the  constraints, 


x  =  Ax  +  br}j  . 


(41) 


By  writing  equation  (40)  in  the  form 


J(tf,ti)  =  J  (x1Qx+n^l,ni+x*'sx+x*Sx)ilt+xT(ti)Sx(t.. )  (42) 

t  . 
i 


where  x(t  .)  is  a  constant,  the  problem  can  now  be  treated  by  the  classical 
calculus  of  vaviat ions . To  remove  the  constraints  imposed  by  equation 
(41),  the  adjoint  vector  of  Lagrange  multipliers  A  is  introduced  and  the 
unconstrained  functional 

G  =  x^Qx+iljPllj+x^Sx+x^Sx+A^  (Ax+brij-x) 

is  constructed,  for  which  the  Euler- Lagrange  equations, 

d/dt  (9G/3X)  ••  3G/9X  =  0, 

where 

VT  _  f  T  ,  T  „  , 

X  -  fx  ,  A  ,  n  j  I  * 

become  x  =  Ax  ♦  bn,  (43) 
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and  the  t ransversal  ity  condition*^  for  free  x(t^)  and  t  Rives 


|i(tf)  =  Sx(tf) 


Substituting  for  i]j  in  equations  (A3)  from  (45)  gives, 


-bP*lbT  x 


These  equations  could  be  solved  for  |j  and  hence  from  equation  (45),  n j  would 
he  known  as  a  function  of  time.  This  is  not  the  form  of  solution  required  for 
feedback  control  where  i)^  is  needed  as  a  function  of  x.  Equations  (48)  show 
that  p  is  a  Linear  function  of  x  and  making  the  substitution 


M  =  Rx 


gives  tlie  matrix  Riccat.i  equation 


R  +  RA  +  ATR  -  RbP-1bTR  -*  Q  =  0 


=  -l>_lblRx 


which  is  now  the  required  form  for  r|j.  The  boundary  condition  for  R  is  found 
from  equations  (47)  and  (49)  to  be 


R(tf.)  a  S 


Equation  (50)  can  thus  be  solved  for  R(t),  t.^t^t  by  allowing  time  to  run 
backwards  and  solving  for  R(t)  where  T  =  t  -t.  This  is  equivalent  to  solving 
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K  +  KA  +a‘k 


-  RI>I>"Vr  1  Q-o 


(53) 


to  fiml  R  as  a  function  of  T  subject  to 


Rt_  «  S 
T=0 


(54) 


Having  solved  for  R,  the  feedback  gain  matrix  K  is  defined  as 


T  -IT 
K  =  P  b  R 


(55) 


and  hence  from  equation  (51) 

i|j  =  -k'x  (56) 

Equation  (56)  implies  that  the  greatest  reduction  in  the  variance  of  the 
perturbed  state  variables  occurs  when  all  the  states  are  fed  back. 

Material  from  the  cited  literature  which  bears  directly  on  liicar  optimal 
control  design  theory  has  been  brought  together  in  the  sections  which  follow 
but  with  some  additional  details  to  emphasise  the  design  techniques  used  in 
this  thesis.  Such  details  concern  the  properties  of  the  Riecati  matrix,  its 
link  with  Liapunov  theory,  its  asymptotic  behaviour,  and  the  reason  why  all 
state  variables  are  fed  back  in  the  optimal  case. 

5.1  Properties  of  Riccati  Matrix  R 

At  this  stage  it  is  worth  examining  the  properties  of  the  Riccati  matrix  K 
a  little  further.  Remembering  that  the  weighting  matrices  Q,  P,  S,  were 
chosen  to  be  symmetric,  then  taking  the  transpose  of  equations  (50)  and 
(52)  shows  that  R*  satisfies  the  same  equation  and  boundary  condition  as  R. 
Thus  R  is  a  symmetric  matrix  and  instead  of  n^  unknown  elements  there  are 


(n-t  1  )n/2 . 
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5.1.1  Relationship  between  R  and  Integral  Performance  Index 

It  is  also  found  that.  R  provides  a  men so re  for  the  minimised  performance 
index.  This  is  readily  shown  by  writing  t  for  in  equation  (40)  and 
substituting  the  optimal  value  -P  ^h*Ux  for  llj  in  equations  (40)  and 
(4t ) ,  giving, 


J(t 


xT(Q+RbP"!hTR)xdT+xT(t  )Sx(t  ) 


(57) 


and 


x  =  (A  -  bP'1brR)x 


(58) 


Differentiating  equation  (57)  with  respect  to  t,  gives 

.T(tf.t)  =  -xT(t)  (Q  •  RhP_1hTlO  x(t)  (50) 

Kxpressing  J(t.  .,t.)  in  the  quadratic  form 

J(tf,t)  =  xT(t)  C(t)  x(t) ,  (60) 

substituting  into  equation  (59),  and  eliminating  x  through  equation 
(58),  yields 

Ct  CA  +  ATC  -  RbP "  ^  hTC  +  (R-C)bP' !bTR  +  Q  =  0  (61) 


where  the  boundary  condi ton 


C(tf)  =  S 


follows  from  equations  (57)  and  (60).  Comparison  of  equations  (61)  and 


(50)  shows  that 


C  =  R 

is  a  solution  of  equation  (61)  and  hence  from  equation  (60) 

x 1 ( t )  K(t)  x ( t )  =  J(tf,t)  (62) 

and  thus  N  provides  a  measure  for  the  minimised  performance  index. 

5.1.2  Link  between  R  and  Liapunov  Stability 

The  Riccati  matrix  R  may  also  he  used  to  check  the  stability  of  the 

/  l  A  A 

optimised  system  through  Liapunov's  second  method.  v  The  scalar 
function  J(tf,t)  can  he  regarded  as  a  Liapunov  function  which  assures 
stability  for  J(lpf)  positive  definite  and  J(t^;t)  negative  definite. 
When  Q,  P,  and  8  are  positive  definite  then  .1(1^1)  must  be  positive 
definite  and  from  equation  (59),  the  additional  requirement  that 
(Q+RhP  h  R)  he  posit  ive  definite  assures  that  ,l(t  .,t)  is  negative 
def in i t  e . 

5.1.3  Asymptotic  behaviour  of  R 

The  asymptotic  behaviour  of  R  as  t  -ku  will  now  be  examined. 
Substituting  equation  (56)  into  equations  (40)  and  (41)  gives 


lf 


J(lf,t.)  =  J  xl(Q+Kl,Ki)xdt>xT(tf.)Sx(tt.) 


(63) 


and 


x  -  (A  -  bK  )x 


(64) 


Thus  if  equation  (64)  represents  a  stable  system,  x(l(.)>0  as  t  ai 


id 


hence 
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»  Constant  us  L  >o. 
f  i  f 

Comparing  this  result  with  equation  (62)  shows  that 

R  -»  constant  as  t^.  -*  <*> 

♦ 

and  hence  R  ■*  0. 

This  result  shows  that  if  equation  (53)  is  used  to  solve  for  R  in  the 
case  of  the  time-invariant,  infinite  upper  limit,  optimal  system,  and  R 
approaches  a  finite  limit  for  large  t  then  the  system  must  be  stable. 
In  this  case  equation  (55)  shows  that  the  optimal  feedback  gain  vector  K 
is  constant. 

5.1. A  Time  invariant,  infinite  upper  limit  performance  integral 

Another  property  of  the  Riccati  matrix  as  determined  from  equation  (50) 
is  that  it  is  independent  of  the  values  of  the  state  variables  and 
hence,  from  equation  (55),  so  are  values  of  the  optimal  feedback  gains, 
K.  This  appears  to  be  a  consequence  of  the  requirement  that  all  state 
variables  be  fed  back.  A  closer  look  at  this  situation  can  be  had  by 
rewriting  equation  (62)  for  the  time-invariant  system  with  infinite 
upper  limit  performance  integral  to  give 

xT(t)R  x(t)  =  J(«,t)  (65) 

In  this  case,  where  the  feedback  gain  vector  K  is  constant,  the  partial 
derivatives  of  J(<*,t)  with  respect  to  each  component  of  K  must  be  zero 
for  a  minimum  value  of  performance  index.  This  means  that  the  partial 
derivatives  of  R  on  the  left  hand  side  of  equation  (65)  with  respect  to 
the  K  components  must  all  be  zero  if  the  minimum  is  to  be  independent  of 
x,  which  evidently  it  is.  Presumably  then,  if  not  all  of  the  state 
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variables  are  fed  hack  or  if  not  aLl  the  feedback  gains  are  regarded  as 
free-  parameter:; ,  not  all  l  he  partial  derivatives  of  the  elements  of  R 
can  simultaneously  be  zero.  Heme  any  sub  optimal  system  based  on 
incomplete  state  variable  feedback  must  consider  initial  values  of  the 
state  variables.  This  result  is  important  because  in  many  practical 
situations  it  is  not  feasible  to  feed  back  all  the  state  variables  and 
there  is  then  a  need  to  examine  sub  optimal  systems.  This  problem  will 
be  discussed  further  in  section  7. 

5.1.5  Alternative  to  Riccati  Matrix  Solution 

As  an  alternative  to  the  Riccati  matrix  equation  solution  for  the  time- 
invariant  infinite  upper  time  limit  optimal  system,  there  is  equation 
(48).  This  equation  together  with  equation  (64)  provides  some  links 
between  the  eigenvalues  of  the  closed  loop  system  matrix  and  the 
elements  of  the  weighting  matrices  P,  Q,  and  the  optimal  feedback  gain 
components  of  K.  Equation  (48)  represents  a  system  of  order  2n  but  it 
is  found  from  symmetry  that  the  characteristic  polynomial  has  only  n+1 
even  powered  terms.  The  characteristic  equation  in  terms  of  the  scalar 
parameter  A  for  this  system  is  given  by  the  determinantal  equation 


\1  -A  b?‘ 1 b1 
n 


~  0 


Q  XI  +A1 

n 


(66) 


where  1^  denotes  the  unit  matrix  of  order  n.  By  virtue  of  the  symmetry 
of  Q  and  bP  b  ,  the  left  hand  side  of  equation  (66)  is  found  to  be  an 
even  function  of  A.  This  fact  is  readily  established  by  substituting  -A 
for  A  in  the  determinant  and  interchanging  rows  and  columns.  Thus  the 
characteristic  polynomial  of  this  optimal  system  is  a  polynomial  in  A2 
of  order  n,  and  the  eigenvalues  occur  in  pairs,  consisting  of  one 
positive  and  one  negative  value.  For  x  to  vanish  as  time  becomes  large 
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the  eigenvalues  with  negative  real  parts  are  assigned  to  x  and  those 
with  positive  real  parts  to  |j.  The  characteristic  polynomial 
corresponding  to  the  x  eigenvalues  is,  from  equation  (64),  the  expansion 
of  the  determinant 

AI  -  A+bK1 

I  n 

and  is  found  to  have  coefficients  which  are  linear  in  the  components  of 
K.  This  polynomial  can  also  be  constituted  from  the  roots  of  equation 
(66)  with  the  negative  real  parts.  Equating  corresponding  coefficients 
in  the  two  polynomials  leads  to  a  linear  set  of  algebraic  equations  for 
determining  values  of  the  components  of  the  optimal  gain  vector  K.  When 
direct  algebraic  factorisation  of  the  characteristic  polynomial  of 
equation  (66)  into  x  and  p  eigenvalue  components  is  possible,  direct 
algebraic  relationships  between  the  matrix  elements  of  A,P,Q,  and 
optimal  gain  components  can  be  obtained.  Having  discovered  that  the 
coefficients  of  the  characteristic  polynomial  in  equation  (64)  are 
linear  functions  of  the  feedback  gain  components  il  is  evident  that  if 
all  the  state  variables  are  fed  back  there  are  sufficient  free 
parameters  to  place  the  system  poles  in  any  desired  positions. 

5.1.6  Selection  of  Weighting  Matrices 

There  now  remains  the  question  of  how  to  choose  the  weighting  matrices 
P,Q,  and  S.  The  ideal  control  system  would  immediately  eliminate 
disturbances  as  they  occurred  but  in  doing  so  would  require  unlimited 
control  power.  To  be  realistic,  finite  errors  due  to  disturbances  must 
be  allowed,  so  that  control  power  demanded  can  be  limited.  The  design 
aim  then  is  to  weight  errors  in  the  state  variables  as  heavily  as 
possible  while  constraining  control  power  to  a  practical  limit  and  at 
the  same  time  ensuring  that  the  assumed  linear  laws  are  not  violated  and 
the  linearised  theory  rendered  invalid.  Dealing  first  with  the 


weighting  matrix  S,  which  because  it  penalises  terminal  errors,  is  used 

in  designing  missile  terminal  control  systems  and  often  leads  to  demands 

for  rapid  changes  in  feedback  gains  when  homing  in  on  a  target.  The 

weighting  matrix  Q  is  used  to  design  systems  which  regulate 

perturbations  in  state  variables  throughout  the  whole  trajectory.  It  is 

this  sort  of  regulation  which  is  sought  in  mid-course  control  of  the 

gliding  projectile  with  consequent  relatively  small  changes  in  feedback 

gain.  Thus  in  the  mid-course  control  system  design  of  this  gliding 

projectile,  S  will  be  taken  as  zero.  A  very  convenient  way  of  dealing 

with  Q  is  to  try,  by  physical  arguments,  to  find  fixed  ratios  between 

its  elements,  and  then  to  take  some  scalar  multiple  of  Q  as  a  variable 

parameter.  Thus  performance  can  be  assessed  against  the  variation  of 

this  one  parameter.  Unfortunately  the  method  cannot  be  applied  to  this 

projectile  because  there  are  absolute  limitations  on  angle  of  attack  a 

and  elevator  control  fj .  A  feature  of  a  high  aspect  ratio  aircraft-like 

wing  is  that  the  variation  in  lift  coefficient  with  angle  of  attack  may 

be  reliably  regarded  as  linear  up  to  the  stall,  but  when  stall  does 

occur,  the  loss  in  lift  is  dramatically  sudden.  Thus  perturbations  in 

angle  of  attack must  be  constrained  to  lie  between  the  steady  value 

a  and  the  stall  value.  If  large  perturbations  in  a.  are  unavoidable  it 

may  be  necessary  to  reduce  a q  and  accept  the  accompanying  reduction  in 

range  in  return  for  improved  accuracy.  In  similar  fashion,  rj j  is 

constrained  to  lie  between  r|Q  and  elevator  deflection  stall  angle.  With 

only  one  control  variable,  P  is  a  scalar,  and  Q  will  have  elements  . 

On  physical  grounds,  reasonable  first  guesses  for  P  and  Q  are  found  by 

choosing  acceptable  maximum  values  for  each  state  variable  perturbation 

x.  ,  and  elevator  angle  perturbation  n,  and  putting  q..  =  l/(x.  x.  ), 
im’  'lm  r  6  nj  lm  jm  ’ 

P  =  1/n2,  Using  these  values,  the  system  equations  can  be  solved  and 
lm 

the  maximum  values  of  the  critical  state  and  control  variables 
determined.  Determination  of  weights  can  then  proceed  by  an  iteraetive 
procedure  which  increases  state  variable  weights  when  full  allowable 


iffpgq  »my V;  ay  wvy^y.  <  y**>  go****  t£« 
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elevator  deflection  is  not  realised  and  vice  versa. 


6.  CONTROL  SYSTEM  DESIGN 

Using  the  physical  properties  of  the  gliding  projectile  given  in  Appendix  I, 

values  of  the  state  variables  V,  a,  8,  q,  h,  X,  were  computed  for  maximum 

range  in  undisturbed  flight  starting  from  an  initial  altitude  of  4  km.  This 

was  accomplished  by  integrating  equations  (6)  to  (11)  using  the  Differential 

(14) 

Analysis  R.  Eplacement  (DARE)  simulation  language  on  the  University  of 

Adelaide's  CDC  6400  computer  and  later  repeated  with  use  of  the  Continuous 
System  Modelling  Program  (CSMP)  on  the  DRCS  IBM  370  machine.  The  set  of 
variables  which  maximise  the  range  in  undisturbed  flight  are  distinguished  by 
the  subscript  0  while  subscript  i  denotes  initial  condition.  Initial  values 
of  Vq  and  0^  are  found  from  equations  (27)  and  (29)  while  aQ  and  n0  are 
obtained  from  equations  (31)  and  (32).  For  the  undisturbed  maximum  range 
trajectory,  equations  (31)  and  (32)  give 

aQ  =  0.17037  rad 
H0  =  -0.55370  rad 

and  for  hQi  =  4000  m,  equations  (27)  and  (29)  give 

V  .  =  172.4  m  s"1 
01 

0  .  =  0.05373  rad 
or 

In  addition,  putting 


allows  equations  (6)  to  (11)  to  be  integrated  to  obtain  the  0  subscripted 
undisturbed  maximum  range  state  variables  corresponding  to  the  constant 
elevator  setting  f|n.  These  results  are  plotted  in  figures  5  to  8,  and  for 


for  projectile  aitr.  no  control  and  control 
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comparison,  some  spot  values  calculated  from  the  approximations  given  l>y 

equations  (26),  (27),  (29),  (30),  are  superimposed  and  show  close  agreement. 

The  0  subscripted  variables  are  the  command  values  of  the  state  variables  that 

the  control  system  must  endeavor  to  achieve  during  disturbed  flight.  Such 

command  values  would  need  to  bo  stored  on-board  in  microprocessor  memory  and 

be  made  available  for  use  by  the  control  system.  The  accurate  generation  of 

command  values  from  the  simple  relationships  given  by  equations  (26),  (27), 

(29),  (30)  presents  the  possibility  of  designing  a  control  system  which  can 

adapt  to  large  errors  in  initial  values  of  state  variables  through  an  on-board 

mini-computer  together  with  memory.  However,  this  paper  is  confined  purely  to 

feedback  control  with  a  preset  store  of  command  values.  The  state  variable 

feedback  scheme  being  consideicd  is  shown  in  figure  9,  where  the  feedback 

gains  Kv  to  are  determined  by  assigning  weights  to  the  five  state  variable 

and  control  variable  perturbations  appearing  in  the  matrices  Q  and  P  of  the 

integral  performance  index  of  equation  (90).  Before  embarking  on  a  design 

procedure  it  is  necessary  to  check  whether  the  desired  linearised  system 

needed  is  both  controllable  and  observable.  When  some  states  are  decoupled 

from  the  input,  the  system  is  not  controllable  and  when  there  is  decoupling 

between  some  of  the  states  and  the  output,  the  system  is  not  observable. 

Tests  for  controllability  and  observability  are  provided  by  the  neat  results 
( l  r> ) 

of  Kalman  .  The  result  for  the  system  given  by  equation  (41)  is  that  for 
controllability,  the  matrix  whose  nth  column  is  given  by  An  ^b  should  he  non¬ 
singular,  where  n  is  the  order  of  the  square  matrix  A.  The  controllability 
matrices  for  the  values  of  A  and  b  given  in  Appendix  1  were  found  to  be  non¬ 
singular. 

6.1  Choice  of  Performance  Index  Weights  and  Determination  of  Feedback 

The  weighting  of  the  state  variables  perturbations  in  the  integral  of  the 

performance  index  will  be  chosen  to  be  of  the  form 
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j  n 

XV  =  2  (x3./xiw)- 

j-1 

and  the  control  weighting  to  he  (H|  /  'ljw)‘: 

Thus  valuer,  of  x.,  n,  which  exceed  x.  n.  ,  will  bo  heavily  penalised.  It 

J  X  J  W  J  w 

will  be  assumed  that  initial  errors  of  20  m  s  1  in  velocity  and  100  m  in 

aJtitude  may  be  encountered.  As  a  first  guess  then,  put 

V,  '  20  m  s"J.  h.  -  100  rn 

Iw  lw 

There  seems  to  be  no  reason  to  limit  the  angular  pitching  rate  q,  hence 
this  term  will  be  omitted  from  the  performance  index.  The  question  of 
maximum  allowable  angle  of  attack  will  now  be  considered.  Maximum  win j 
lift  is  dependent  on  Mach  number,  aspect  ratio,  and  wing  cross-section. 
Design  handbooks  such  as  references  7,  S  give  values  of  maximum  lift 
coefficient  for  various  Mach  numbers  and  aspect  ratios.  Tor  a  wing 

of  aspect  ratio  7  at  a  Mach  number  of  0.5,  which  corresponds  roughly  to 
V  .  172.4  rn  s  J  at  altitude  h  .  ~  4  km.  a  CruAV  of  about  1.4,  based  on 
wing  area,  should  be  achievable.  The  maximum  range  flight  at  angle  of 
attack  -  0.17037  radian  requires  to  be  1.16.  This  leaves  a  margin  of' 
about  0.035  radian  or  2  degrees  for  control  purposes.  The  constant 
elevator  angle  deflection  for  maximum  range  flight  is  -0.5537  radian  so 
that  if  a  maximum  deflection  of  about  0.7  radian  or  40  degrees  is  permitted 
a  reasonable  choice  for  q,  is  0.1  radian.  An  initial  guess  for  an 

attitude  angle  weight,  0  ,  can  be  obtained  by  considering  the  flight  paths 

|W 

of  figure  2.  For  100  m  error  in  altitude  at  the  end  of  midcourse  flight 
the  corresponding  mean  flight  path  error  is  about  1/4  degree  and  the  change 
in  line  of  sight  angle  to  target  is  about  1  degree.  Allowing  4  degrees  of 
total  seeker  beam  angle  to  accommodate  final  midcourse  altitude  errors 
requires  a  mean  midcourse  flight  path  error  of  i/2  degree.  Thus  as  a  f;u  >1 

guess,  0  will  be  taken  to  be  1/2  degree.  Summarising,  the  weights  on  the 

lw 


stale  and  control  vocable  perturbations  wore  first  taken  to  bo 


0.035  r^d.ian 


0  -  0.0087  radian 

jv 

h .  -  100  in 

!u 

n,  0.1  radian. 

!  w 


l'ho  weighting  matrices  0  and  I'  appearing  in  equations  ;-!0)  and  (5J)  at 
t  ints  taken  as 


l/Vlw 


1  /c; .  ' 
Iw 


and  1'  l/n,;> . 

'  Mv 


The  control  system  resulting  from  this  choir*  of  0  and  P  will  be  toferrod 

to  as  System  I.  substituting  these  values  for  Q  and  V  into  the  KiccaU 

matrix  equation  (55),  together  with  values  of  a  and  b  obtained  iron 

equations  v  >-l )  to  (IS),  allowed  R  to  be  determined  and  hence  the  optimal 

feedback  gain  vector  K  from  equation  (55).  Because  the  Riccati  matrix  in 

this  case  is  of  order  >x;>  and  is  symmetrical  there  are  fifteen  unknown 

elements  r  .  to  bo  solved  for  from  the  set  of  simultaneous  equations 
fl 


r  a,  .  i  .•  i .  .a.  .  i  q . 
i|  kj  kj  Kl  ‘  l  j 


' '’Iw'  ‘Vo  j 1  i  j '  b2b-» ( ' 12r -t  j  ' rMr2j 1 ’ b7  L-l r  i  j  ! 
where  for  j  -  I ,  i  I,')  ,1.1, 5 
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j  ^ ,  1  "■  ^ » 4 ,5,6, 
j-3,  »  =  3,4,5, 
j -4,  i  =  4,5 
j=5,  i  =  5. 


and  tin*  iii'c  elements  of  Q.  Solutions  tor  the  r .  .  were  computed  from 
equations  (68)  using  both  the  DARK  and  CSMP  integration  routines. 
Components  of  the  time  varying  otpimal  gain  vector  K  were  then  found  from 
equation  (55)  and  equation  (56)  gave  the  elevator  deflection  angle 

increment  i^,  which  when  added  to  q^,  gives  the  value  for  i)  to  he 

substituted  into  the  trajectory  equations  (6)  to  (ll).  These  equations 
were  integrated  for  an  initial  velocity  20  in  s  *  greater  than  t ho  initial 
optimal  value  of  V  .  -  172.4  m  s  V  The  resulting  trajectories  for  V,  0, 
h,  are  shown  in  figures  5  t.o  8  where  the  optimal  and  no-feedback  control 
trajectories  are  shown  for  comparison.  These  figures  clearly  reveal  the 
low-damped  oscillatory  phugoid  motion.  From  figures  7  and  8,  the  optimal 
flight  path  has  an  altitude  of  1  km  when  the  ground  range  is  25  Um  and  the 
time  taken  for  the  projectile  to  reach  this  position  is  158  s .  The  Very 
marked  effect  that  velocity  perturbation:;  have  on  attitude  is  indicated  by 
figure  6,  which  shows  that  the  initial  velocity  perturbation  of  20  m  s  * 
induces  an  attitude  oscillation  of  about  7  degrees  amplitude  initially  in 
uncontrolled  Might,  decaying  to  1.7  degrees  at  t  =  158  s.  In  controlled 
flight  the  maximum  amplitude  is  2.4  degrees,  decaying  to  0.7  degree. 
Figure  7  shows  how  the  controlled  flight  improves  the  altitude  variation 

with  time  in  that  the  rate  of  descent  is  closer  t.o  the  optimum  than  the 

uncontrolled  flight.  However,  it  is  the  positional  accuracy  which  is 
important,  and  figure  8  reveals  that  the  controlled  I  light  shows  no 
improvement  over  uncontrolled  flight  at  a  ground  range  of  25  km.  The 
altitude  error  at  25  km  ground  range  is  about  300  m  which  would  require  a 
seeker  beam  angle  of  6  degrees  t.o  accommodate,  and  to  which  must  he  added 
another  1.4  degrees  for  the  attitude  error.  At  this  point,  the  question  to 
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be  investigated  is  whether  this  seeker  beam  angle  requirement  of  7.4 
degrees  to  accommodate  an  initial  velocity  perturbation  of  20  m  s  *  (-an  be 
reduced  by  improving  the  feedback  control. 

To  this  end,  further  changes  in  weighting  values  were  made.  Values  of  a 
and  r)jw  were  left  unchanged  because,  in  the  trajectories  just  described, 
values  of  maximum  angle  of  attack  and  elevator  control  deflection  were 
considered  to  be  marginally  acceptable.  The  procedure  was  then  to  diminish 
the  weights  ,  0^,  h^,  by  a  factor  of  two  successively  and  to  calculate 
trajectories  for  various  combinations  of  these  diminished  weights  and  to 
note  changes  from  optimal  terminal  conditions  at  t  =  158  s  in  response  to 
an  initial  velocity  perturbation  of  20  m  s  Because  a  check  had  shown 

that  there  was  not  a  significant  difference  between  trajectories  calculated 
with  time-varying  optimal  feedback  gains  and  those  calculated  for  constant 
feed  back  gains  based  on  values  of  Riccati  elements  corresponding  to 
T“!58  s  such  constant  gains  were  used  in  the  trajectory  selection  procedure 
to  reduce  excessive  computation.  To  assist  in  the  selection  process  it  is 
possible  to  compute  and  compare  values  of  the  performance  index  given  by 
equation  (40),  for  the  various  feedback  gain  vectors.  This  is  achieved  by 
solving  equation  (68)  for  the  r„  with  the  r)^2  term  omitted  and  by 
substituting  a..-b.K.  for  a  .  and  q . .fK.K./n,2  for  q...  The  underlying 
philosophy  behind  this  approach  is  that  a  tightening  of  the  control  over 
each  state  will  bring  each  state  close  to  the  required  optimal  value  at 
corresponding  times  and  this  should  then  also  lead  to  acceptable  position 
accuracy.  However,  attempts  to  exercise  tight  control  over  attitude  and 
velocity  simultaneously,  revealed  some  adverse  interaction  arid  a  closer 
examination  of  the  physical  reasons  was  demanded.  This  adverse  interaction 
st. nis  from  trying  to  regulate  velocity  through  elevator  controls  without 
significantly  changing  the  projectile's  attitude  or  flight  path  direction. 
To  regulate  a  velocity  in  excess  of  the  desired  value,  the  elevators  would 
be  deflected  so  as  to  increase  the  angle  of  attack  and  hence  increase  the 


induced  drag.  Thus,  during  a  ptiugoid  velocity  cycle,  the  higher  velocity 
purl  would  lie  Mown  at  higher  than  optimum  angle  ol  attack  and  the  lower 


velocity  part  at  lower  angle  of  attack.  In  this  way  the  mean  drag  force  is 
increased,  as  is  required,  but  so  also  is  the  mean  lift  force  and  it  can  be 
shown  that  the  undamped  amplitude  of  the  phugoid  flight  path  oscillation  is 
approximately  proportional  to  the  square  root  of  the  lift  coefficient  C^. 
This  result  is  readily  derived  by  solving  the  open  loop  part  of  the  reduced 


system  equations  (11.35)  for  with  z^O  and  z^=Vj.  at  zero  initial  time. 
In  this  notation,  z^  is  velocity  perturbation  and  is  the  flight 
perturbation  angle  iji^  .  It  is  found  that  the  undamped  amplitude  of  ijjj  for 


the  above  initial  conditions  is  proporLional  to 
(35), 


'a21/uV 


From  equation 


a21  3  -(prefS/2m)CI (a0)-(g/v2)cos(00-a0)  (69) 

when  the  contribution  of  elevator  lift  is  ignored.  Under  optimal 
conditions,  lift  is  approximately  equal  to  weight  and  hence  the  magnitude 
of  the  two  terms  in  equation  (69)  are  roughly  equal  so  that 


a21  ”  -2(PrefS/'2m)c:j  (w0) • 


Equation  (11.32)  shows  that.  ui.  is  proportional  to  (C,  (a„))  '  and  thus 

ip  L  0 


-a^j/u)^  is  approximately  proportional  to 


<w)l/2- 


Hence  attempts  to  control  velocity  in  this  way  tend  to  accentuate  the 
phugoid  oscillation  initially  by  increasing  the  frequency  and  undamped 
amplitude  of  the  flight  path  angle  although  damping  is  increased.  This 
initial  effect  tends  to  increase  the  velocity  direction  error  and  although 
the  velocity  magnitude  error  may  be  reduced,  a  flight  path  positional  error 
is  still  ma i nta  ined . 


Nov;  consider  what  happens  when  the  elevators  act  in  the  opposite  manner  so 
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that  they  reduce  angle  of  attack  in  response  to  an  increase  in  velocity. 
The  tendency  then  is  to  reduce  tiie  phugoid  frequency  and  undamped  amplitude 
and  although  the  velocity  magnitude  is  not  so  well  regulated,  the  flight 
path  positional  accuracy  tends  to  improve.  However,  if  this  form  of 
control  is  overdone  in  such  a  way  that  angle  of  attack  is  reduced  too  much 
in  response  to  increased  velocity,  an  instability  will  be  reached  when  the 
glider  will  begin  to  nose  dive  at  ever  increasing  velocity. 

This  interaction  between  velocity  and  attitude  control  of  a  glider 
highl Lghts  the  difference  between  flight  path  control  of  a  powered  aircraft 
and  a  glider.  On  the  powered  aircraft,  velocity  regulation  would  be  sought 
through  feedback  to  engine  throttle  as  well  as  to  elevator  controls. 

Returning  now  to  the  search  for  more  accurate  trajectories  by  diminishing 

V,  ,  (-)  ,  h ,  ,  it  was  found,  in  accordance  with  the  above  considerations, 
lw’  iw’  lw 

that  positional  accuracy  was  not  improved  by  reducing  vlw>  but  that 
improvement  was  gained  by  reducing  0^  from  1/2  degree  to  1/4  degree  and  h^ 
from  100  m  to  12. 5  m.  For  the  set  of  weights 

V,  =  20  m  s"1 
lw 

a,  -  0.03T>  radian 

lw 

0  =  0.00436  radian 

iw 

h ,  =  12.5  m 

lw 

n,  =0.1  radian 

lw 

equations  (68)  and  (55)  were  again  used  to  obtain  an  optimal  time  varying 
gain  vector  K.  The  components  of  K  are  plotted  in  figure  10.  This  figure 
shows  the  gain  components  to  be  virtually  constant  over  most  of  the  flight 
time  and  consequently  trajectories  were  computed  for  these  gains  and  for 
gains  held  constant  at  their  initial  values.  These  trajectories  were 
sufficiently  close  to  each  other  to  indicate  that  constant  feed  back  gains 
based  on  a  constant  initial-valued  plant  matrix  A  and  control  vector  b  were 
adequate  for  design  purposes.  This  is  a  useful  attribute  of  the  optimal 
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control  approach  in  that  it  does  indicate  when  a  time-invariant  system 
representation  is  a  reasonable  approximation  for  the  determination  of  state 
feedback  gain  and  thus  when  the  extra  complication  of  time  varying  gains 
can  be  dispensed  with. 

The  values  of  feedback  gains  used  in  the  constant  gain  trajectory 
computations  are:- 

Ky  =  -2. 034  x  10  2  radian  m  *s 
Ka  =  6.335 
=  -26.469 

0 

K  =  -1.867  s 

q 

K,  =  -7.566  x  10  3  radian  m  * 
h 

The  control  system  using  these  values  of  feedback  gain  will  be  referred  to 
as  System  II.  Trajectories  computed  for  this  set  of  gains  are  shown  in 
figures  1 1  to  16  for  an  initial  velocity  perturbation  of  20  m  s  1 .  From 
figure  14  the  altitude  error  at  25  km  ground  range  is  200  m  and  from  figure 
12  the  attitude  error  at  t  =  158  s  is  0.45  degree.  Thus  the  seeker  beam 
angle  needed  to  accommodate  this  error  is  4,9  degree  or  0.245  degree  per 
n  s  *  of  initial  velocity  error.  In  figure  12  the  oscillation  in  0Q  is  due 
to  the  magnified  scale  showing  up  round  off  error  in  the  value  used  for  Vq. 
and  stresses  the  sensitivity  of  attitude  angle  to  velocity  perturbations. 
Figure  13  shows  that  the  controlled  rate  of  descent  is  indistinguishable 
from  the  optimal  value,  and  the  altitude  error  shown  in  figure  14  arises 
from  the  positive  velocity  increment  shown  in  figure  11.  In  achieving  this 
amount  of  control,  figure  15  indicates  a  maximum  angle  of  attack  variation 
from  optimum  of  0.0322  radian  or  1.84  degree  and  from  figure  16  the  maximum 
elevator  deflection  angle  is  0.515  radian.  Hence,  in  responding  to  an 
initial  velocity  error  of  20  hi  s  '  the  angle  of  attack  is  contained  to 
within  2  degrees  of  the  optimum  value,  as  required,  and  the  elevator 
deflection  angle  remains  less  than  0.7  radian.  Another  trajectory  which 
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was  computed  for  an  initial  altitude  error  of  100  in  showed  that  this  error 
had  been  halved  by  the  end  of  midcourse  flight. 

To  demonstrate  the  stability  of  this  control  system  to  in-flight 
disturbances,  a  series  of  random  wind  gusts  was  introduced  along  the 
trajectory  having  an  initial  velocity  error  of  20  m  s  *  Horizontal  and 
vertical  components  of  these  gusts  are  shown  in  figures  17  and  18,  and 
figure  19  shows  them  superimposed  on  projectile  velocity.  The  gusts  were 
of  duration  5  s  and  had  random  directions  and  random  magnitudes  in  the 
range  0  to  10  ms  *.  Responses  to  these  gusts  and  to  the  initial  velocity 
error  of  20  m  s  of  attitude,  altitude,  angle  of  attack,  and  elevator 
deflection  angle,  are  shown  in  figures  19  to  24 ,  and  are  all  stable. 

6.2  Integrity 

There  is  now  the  question  of  the  integrity  of  this  control  system  11  and 
the  effect  of  the  single  failure  of  each  state  variable  feedback  will  be 
examined.  In  a  time- varying  linear  system  one  could  compute  the 
performance  index  by  modifying  the  Riccati  equation  in  the  manner 
previously  indicated  in  the  search  for  a  suitable  weighting  matrix  Q.  Such 
computations  readily  reveal  a  non-converging  asymptotic  behaviour  when  It 
exists  and  can  thus  be  used  to  indicate  stability  or  the  lack  of  it. 
However  in  this  instance  where  the  system  is  well  approximated  by  a  time- 
invariant  system  it  is  easy  to  compute  eigenvalues  of  oloseu-loop  matrices 
formed  from  the  difference  between  the  ini  t  iai-valu,  plant  matrix  and  the 
product  of  the  initial-valued  control  vector  and  a  constant  feedback  gain 
row  vector  which  has  each  gain  component  made  zero  successively. 
Eigenvalues  of  the  closed-loop  system  IT  having  feedback  failure  of  one 
state  at  a  time  are  given  in  Table  1. 

Table  1  shows  that  attitude  feedback  failure  is  completely  destabilising. 
The  main  effect  of  each  other  failure  is  to  greatly  reduce  short  period 
damping.  Trajectory  computations  for  each  failure,  other  than  attitude, 
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TABLE  1.  SYSTEM  II  POLES  FOR  EACH  FEEDBACK  FAILURE 


Feedback 

Angle 

Attitude 

Fa i lure 

None 

Velocity 

of 

Attack 

-10.62  ±  12 . 23i 

-0.1102  ±  16 . 22i 

-0.2008  ±  18.29i 

-6.1629 

System 

-0.3174 

-0.3148 

-0.2250 

5.618 

Poles 

-0.06582 

-0.07605 

-0.08051 

0.3062 

-1 

s 

-0.006205 

-0.005409 

-0.005689 

-0.3690 

-0.008212 

Feedback 

Pitch 

Altitude 

All 

Failure 

Rate 

-0.1135  ±  16.22i 

-0.1135  ±  16 . 23i 

-0.3456  ±  6 . 045i 

System 

-0.3178 

-0.3774 

-0.01041  +  0 . 08546 i 

Poles 

-0.06555 

-0.01051 

-0.0008673 

“1 

s 

-0 . 006203 

.. 

-0.001628 

confirmed  that  flight  stability  was  maintained,  except  for  pitch  rate 
feedback  failure  in  the  face  of  velocity  perturbations.  For  pitch  rate 
feedback  failure,  angle  of  attack  perturbations  reached  such  large 
amplitude  that  for  negative  perturbations  the  lift  became  significantly 
less  than  the  weight,  causing  the  projectile  to  plunge  downwards  and 
further  increase  the  velocity  perturbation  to  values  well  outside  the 
validity  of  the  linear  theory.  Under  these  conditions,  computations  of  the 
response  to  a  10  n  s  '  initial  velocity  perturbation,  showed  that  the 
projectile  plunged  downwards  and  made  a  180  degree  backward  turn. 
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6.3  Implement a lion 

This  feedback  system  requires  continuous  measurement  of  the  live  state 
variables,  velocity,  angle  of  attack,  attitude,  pitching  rate,  and 

altitude.  One  could  use  a  pitot  tube  for  velocity,  a  yawmeter  for  angle  of 
attack,  a  position  gyroscope  for  attitude,  and  an  altimeter  for  altitude. 
For  pitching  rate  measurement  there  are  rate  gyroscopes,  or  an  array  of. 
accelerometers,  or  coriolis  devices  cuing  deflection  of  small  gas  jets  or 
vibrating  simply  supported  beams.  A  practical  difficulty  arises  from  the 
accuracy  required  for  the  initial  setting  value  of  a  position  gyroscope 
used  to  measure  attitude.  As  mentioned  previously,  a  1/4  degree  error  in 
flight  path  direction  leads  to  a  100  m  error  in  terminal  altitude  which  in 
turn  requires  2  degrees  of  seeker  beam  to  see  the  target.  11  the 
projectile  is  launched  from  a  moving  platform  such  as  an  aircraft  then  an 
accuracy  of  a  fraction  of  a  degree  in  the  initial  gyroscope  axis  sotting 
would  not  be  achievable  by  simple  means  and  attitude  would  not  be  available 
for  mid  course  trajectory  correction.  When  states  are  nut.  available  for 
measurement,  the  sophisticated  solution  is  to  use  an  on-board  computer 
simulation  of  the  plant,  behaviour  to  estimate  the  missing  states  from  the 
measured  ones.  The  way  in  which  this  is  usually  carried  out  is  described 
with  reference  to  figure  2T> .  This  figure  shows  the  output  from  the  part,  oi 
the  system  under  control  to  be  the  state  vector  x.  When  not  all  the 
components  of  x  are  available  for  measurement,  the  remaining  ones  may  he 
represented  by  Fx  where  F  is  a  matrix  which,  for  example,  would  be  a 
unitary  matrix  if  all  of  x  was  available.  Thus  the  signal  which  passes  on 
into  the  estimating  part  of  the  scheme  consists  of  the  measurable  parts  of 
x  contaminated  by  the  measurement  noise  v.  From  this  signal  an  estimate  ot 
the  state  vector,  denoted  by  x,  is  computed  and  led  back  to  the  system 
contro'ler.  Clearly  then,  the  function  of  the  estimation  process  i  to 
make  x  approximate  x  as  closely  as  possible.  From  figure  2!i  it  is  easily 


shown  that 
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x  -  x  =  (A-bK^F) (x-x)  -bK^.v+w  (70) 

This  equation  shows  that  estimation  performance  depends  on  the  choice  of 
the  gain  vector  K^.  For  a  time-invariant  system,  classical  design 
techniques  can  be  used  to  select  K^.  such  that  the  estimation  process  is 
stable  and  such  that  a  quirk  response  minimises  the  difference  between  x 
and  x.  If  the  components  of  the  noise  vectors  v  and  w  are  assumed  to  be 
white  and  Gaussian  with  zero  me*"  values  then  a  sophisticated  method  of 
determining  K^.  is  available  through  use  of  Kalman  filter  theory, 
descriptions  of  which  are  to  be  found  in  references  11  and  16.  The  result 
of  this  theory  is  that  the  best  estimate  of  x  is  given  by  x  when 

k’J  =  P-1FR  (71) 

and  R  is  the  solution  of  the  Riccati  matrix  equation  (50),  integrated 
forward  in  time,  with  the  elements  of  the  matrices  Q  and  R  now  representing 
the  power  spectral  densities  ot  the  components  of  the  system  noise  vector  w 
and  the  measurement  noise  vector  v  respectively. 

Another  approach  to  design  when  some  states  are  not  available  and  when  the 
system  is  time-invariant  is  to  investigate  the  use  of  compensating  networks 
to  achieve  the  closed  loop  transfer  functions  derived  by  optimal  control 
theory  for  complete  state  feedback.  Techniques  for  designing  such 
compensating  networks  are  discussed  in  references  10  and  17.  However,  the 
approach  that  will  be  considered  in  more  detail  in  the  next  section  is  that 
of  approximating  the  system  to  be  controlled  by  one  of  reduced  order. 

7.  SUB  OPTIMAL  CONTROL 

In  this  section  the  possibility  of  sub  optimal  control  is  pursued  because  the 
lower  performance  from  partial  state  feedback  is  compensated  by  design 
simplicity.  The  first  method  to  be  considered  is  that  of  representing  the 
system  by  a  reduced  order  mathematical  model  and  using  this  model  as  the  basis 
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for  applying  linear  optimal  control  theory.  Provided  that  this  model  is  a 
close  approximation,  we  might  expect  the  design  to  he  nearly  optimal  in  the 
sense  that  the  optimality  will  be  almost  independent  of  initial  conditions. 

7.1  Reduced  Order  Approximation 

Reduced  order  approximations  are  likely  to  succeed  when  there  is  wide 
separation  between  system  poles  and  the  concept  of  fast  and  slow  variables 
can  be  invoked.  in  dynamical  problems  this  situation  leads  to  the  use  of 
singular  perturbation  techniques  which  can  be  particularly  useful  in 
solving  non  linear  equations.  A  discussion  on  the  application  of  singular 
perturbation  techniques  to  solution  of  aircraft  performance  problems  is 
given  in  reference  18  where  conditions  for  successful  application  arc 
investigated.  Typically,  for  conventional  aircraft  configurations,  the 
phugoid  and  short  period  open  loop  modes  can  be  separated,  and  the 
projectile  considered  here  is  no  exception  as  is  shown  in  Appendix  11. 
However,  in  using  this  separation  property  to  design  closed  loop  control 
systems  it  is  necessary  to  check  the  effect  of  loop  closure  on  pole- 
positions.  At  the  end  of  Appendix  II  it.  is  shown  that,  the  phugoid  motion 
can  he  approximated  by  the  second  order  system 


zi  =  aiiVai:iWi 

K  =  *a2iVa2.:;VVi 


(72) 


where  z.j  = 


-  0-tKj  =  the  increment,  in  flight  path  angle, 

dl  =  a!2+al3 
d2  =  _(a22+a23)' 


the  elements  a.,  being  those  corresponding  to  the  elements  of  the  complete 
plapt  matrix  A  based  on  initial  flight,  conditions.  Hy  putting 
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iiiul  V  =  l/«.2 

lm 

in  equation  (66),  the  characteristic  equation  for  an  optimally  controlled 
system  based  on  equations  (72)  and  their  adjoint  set,  becomes 
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thus 

A2  -  (A) 

where 

'S 

,  A2  are 

(79) 


Hence  from  equations  (75)  and  (76) 


2\ j  -  (-c^c,1'2)1'2  .(-cflcV2)1'2 


(«0) 


and 


1/2. 1/2  ,  1/2. 

2A2  =  -(-Cj-20  )  -(-Cj+2c3  ) 


1/2 


(81) 


Equations  (75)  to  (78)  provide  the  links  between  the  poles  ,  A.,  of  the 

optimal  system  and  the  assigned  weighting  ratios  |i,  V. 


Optimal  feedback  gains  k(  and  can  be  related  to  the  poles  and 

weighting  ratios  p,v  by  first  substituting 


di =  ’kivrk2+i 


(82) 


in  equations  (72)  so  that  the  characteristic  equation  becomes 


A2  -  (aj  1-a33-d1kl*d2k,)A 

+  ^a2lal  3”ai  ia23^+^l  ^la23+t*2a13^  ”k2^2al  l+<*  1**2 1 


(83) 


Equating  coefficients  in  equations  (79)  and  (83)  and  solving  for  k.  and  k., 


kI  “  ‘'u(l  l*'2l  +  ' 2a\r ’■,21(<,r,23"'<?an) 

^(dla2?Vu),dl(VllJV2iy 

-  (Aj+A^Kd  |tt, ,  M.,a  ( , )  H^A^ 

d2 1  d  1 a  2  3+ d  2°  1 3  V  d  2  a  1  l +  <  V  2 


k2  =  al3(dl;’21+d2all)'a23(<,ia23+,V!3) 

~  i ; 

-  (Aj+A^)  (^ja2.3+<i2a  1 3^”dl\^'>  ) 

*^2  ^ d  1  *'  2  3^  d  2*1 1 3  ^  * d  1  ^  d  2  *  11  ^ d  1 a  2 1  ^ 

Thus  for  the  second  order  system  Riven  by  equations  (72),  plots  >1  optimal 
feedback  Rains  k  and  k0  and  of  weight  .ratios  |.i  and  v  can  be  drawn  in  the  A 
plane  from  the  relationships  given  by  conations  (77'),  (78),  (<U)i,  ioi). 

(HA)  and  (85).  Such  plots  are  shown  in  figure  2 6  lor  tin*  case  when  the' 
characteristic  equation  for  the  closed  Loop  sj,  item  has  a  pair  of  complex 
conjugate  roots.  From  equations  (84)  and  (85)  it  is  seen  that  when  Aj,  A,, 
are  complex  conjugates  of  each  other,  the  looii.  of  and  k,  in  the  A  plane 
are  concentric  circular  arcs.  Putting  A-O+iU),  the  centre  of  the  kj  circles 


at  a  •-  ( ‘1  •>** f ^  | a •>  i ) / ‘I .) 


and  of  the  k.>  circles 


at  o  =  ■(<,i;*'>3f,,2a13)/,,l 


Figure  26  shows  curves  of  constant  values  oi  |.i  along  which  V  varies.  The 
largest  value  shown  tor  p  i n/450  radian  in  s  which  would  correspond  i  ■> 
ii  ~2  degrees  and  V  =5  m  ■-  .  From  equal  ions  (7  7)  and  (78)  it  can  In*  seen 


-  47 


tli.it.  tilt-  roots  of  the  open  loop  ( tin  inc  tor  i  si  i  c  equation  occur  when  v-~|i~0, 
and  at  that  point  on  figure  26  the  values  of  a  and  u>  are  in  close  agreement 
to  those  given  for  the  complete  open  loop  system  in  Appendix  I.  As 
velocity  weight  is  Lightened,  or,  that  is,  V  is  decreased  and  hence  p 
increased,  figure  26  reveals  that  the  phugoid  becomes  more  oscillatory. 
This  is  because  the  system  is  trying  to  regulate  velocity  through  the 
effect  of  angle  of  attack  on  drag.  As  attitude  weight  is  tightened,  that 
is  as  is  made  smaller,  v  increases  and  the  system  response  eventually 
becomes  aperiodic.  In  fact  as  one  root,  of  the  characteristic  equation 
tends  towards  the  centre  of  the  circles  and  the  other  towards  -d^v. 
Relaxing  the  velocity  weight  completely  is  represented  by  the  p=0  curve  in 
figure  26  and  along  this  curve  k^  is  positive,  showing  that  when  velocity 
is  free  to  change,  angle  of  attack  is  decreased  when  velocity  perturbations 
are  positive.  Thus  the  physical  features  of  the  phugoid  response  are  borne 
out  in  figure  26. 

To  relate  the  feedback  gains  ,  k,,  of  the  reduced  system  to  feedback  gains 
in  the  full  system,  the  assumption  is  made  that  is  a  fast  variable  and 
will  respond  instantly  to  elevator  deflection  so  that  the  steady  state 
relationship  of  equation  (24)  can  be  used  to  give 

n,  =  -WjOcyaoO/oc^/bn) 


-  ’  Mi  j  /  -i .  from  Appendix  1 


Thus,  put  l  iug  e  j  -k^'itj 


*  -Vi  “WV 


leads  to  Kj  -  13k, /4,  —  ")3k^.4  K ,,  -  1 3  k  ^  /  o  , 


where  K  ,  ;. ,t .  K,  are  opt  i 


^ ;  *  .  > 


per  il  i  i  j  on  t  oedbuk  of 


V,.  t<  r)  t.c.'  tf  e  1  e>.  <t  ,oi 


48 


Til  order  to  restrict  the  total  elevator  deflection  angle  to  less  than  0.7 

-3  - 1 

radian,  the  magnitude  of  k^  must  not.  exceed  2.25x10  '  radian  ns  s  to 

accommodate  a  maximum  initial  velocity  error  of  20  m  s  *.  Thus  in  figure 

26  the  choice  of  optimal  feedback  gains  must  be  restricted  to  values  lying 

-3  - 1 

in  the  region  to  the  right  of  the  k^  -  -2.25x10  radian  m  s  circle,  which 

is  seen  to  be  approximately  given  by  (J=rr/900  radian  m  1s.  This  value  of  p 

corresponds  to  V,  -10  m  s  ^  for  «,  =2  degrees. 
r  lm  lm  n 

For  aperiodic  response  wlien  the  two  roots  of  the  characteristic  equation  of 

the  reduced  system  lie  along  the  negative  real  axis  the  corresponding 

optimal  feedback  gains  and  weigh*-  ratios  are  plotted  against  the  root 

closest  to  the  origin  and  the  plots  are  shown  in  figures  27  and  28.  As  v 

becomes  large  it  is  found  that  the  root  farthest  from  the  origin  on  the 

negative  real  axis  behaves  as  -d^V,  k^  approaches 

-3  -1 

-a9 j/dn=l . 72x10  '  radian  m  s  and  k^  approaches  v.  Figure  28  shows  that 
for  the  aperiodic,  response,  k^  is  very  close  to  its  asymptotic  value  V. 
Thus  for  large  values  of  V  when  one  of  the  phugoid  poles  moves  far  away 
from  the  origin,  the  assumption  that  the  phugoid  poles  and  the  short  period 
poles  are  well  separated,  may  be  violated  and  a  check  on  the  stability  of 
the  complete  closed  loop  system  with  feedback  gains  and  k.;  should  be 
made.  This  check  showed  that  for  values  of  v  between  2  and  4,  the  feedback 
gains  kj  and  caused  undamping  of  the  short  period  mode  thus  rendering 
the  complete  system  unstable.  The  stability  boundary  corresponding  to  this 
nndnmping  of  the  short  period  mode  is  shown  in  figures  27  and  28. 

There  now  remains  t..e  problem  of  selecting  a  pair  of  feedback  gains  and 
,  which  reduces  to  choosing  relative  weights  |j  and  v.  In  accordance  with 
previous  considerations  that  or  should  be  restricted  to  about  2  degrees  and 
that  the  average  error  in  flight  path  direction  should  be  no  more  than  1/4 
degree,  the  values 


a 


In 


-  n/’J 0  radian 


roa 


-  49  - 


and  ill.  —  K/720  radian,  are  taken. 

1  in  *  ’ 

g  i  v  i  1 1  k  V  ”  8  . 

Because  Vjn(  -  20  in  s  *  was  found  to  give  the  best  results  for  the  complete 
optimal  control  system,  this  value  will  be  selected  again,  giving 

p  =  n/1800  radian  m  ^s. 

But,  as  shown  in  figures  27  and  28,  the  complete  system  is  unstable  for 
v=8,  so  some  other  criteria  must  be  used  for  feedback  gain  selection.  The 
approach  used  here  is  to  search  the  sLable  regions  of  figures  27  and  28  for 
values  of  the  wei.«'r,<  ratios  p  and  v  which  minimise  an  integral  performance 
index  for  the  reduced  system.  In  the  general  notation,  the  quadratic 
integral  performance  index  based  on  the  reduced  system  given  by  equations 
(72),  is 

Jft)  -  z'kz 
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gains,  we  can  put 


U]  =  l~^2Z2  ’n  eVatlon  (86) 

and  letting  t>0  in  the  integral  performance  index,  the  Riccati  equation  for 


R  reduces  to  a  set  of  linear  equations  in  the  elements  r  „  ,  given  by 
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In  solving  these  equations,  q  ^  and  q^  ate  taken  as  the  inverse  squares  of 


the  maximum  allowable  values  of  and  i|i^ ,  which  are 


Vlm  =  20  m  s  ,  »jf lm  =  1/^  degree , 


but  k  and  correspond  to  values  of  p,  V  selected  in  the  region  being 
searched.  When,  as  in  previous  considerations,  only  initial  values  of 
velocity  perturbations  are  considered,  then  from  equation  (86) 
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t+0 
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and  hence  can  be  used  as  a  measure  of  performance  index  .1 .  Such  values 
of  were  calculated  along  the  curves  of  constant  p  shown  in  figures  26, 
27  and  28  except  for  p=n/450  radian  m  's  which  requires  too  much  control 
deflection.  The  lowest  values  of  r  were  found  to  occur  along  p=0, 
decreasing  monotonica 1 ly  from  v=0  to  a  minimum  between  v= 2  and  v~U .  Thus 
values  of  k^  and  are  sought  which  correspond  to  a  point  on  p-0  in  figure 
27,  lying  to  the  left  of  the  stability  boundary  at  a  distance  which  ensures 
adequate  damping  of  the  short  period  inode.  Trajectories  were  computed 
using  the  optimal  feedback  gains  corresponding  to  v~ 1 ,  1.4,  1.6,  and  p-0. 
As  expected,  the  greater  the  value  of  v  the  less  the  altitude  error  but  the 
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lower  the  short  period  damping.  The  trajectory  for  v=l  showed  no 
appreciable  altitude  correction,  whereas  trajectories  for  v=1.4  and  1.6 
both  showed  about  50  in  correction.  Trajectories  for  v=1.4,  which  have  more 
short  period  mode  damping  than  those  for  V-1.6,  are  shown  in  figures  29,  TO 
a  nd  3 1 . 

At  this  stage  it  is  worth  investigating  how  a  similar  approach,  through  use 

of  the  reduced  order  system,  could  be  used  to  achieve  a  tighter  regulation 

of  velocity  at  the  expense  of  flight  path  direction.  This  means  choosing 

Vjm  to  be  very  small  and  tji  to  he  large.  In  terms  of  relative  weights,  |J 

is  to  be  as  large  as  possible  with  v  small.  From  figure  26,  the  largest 

value  of  |.i  which  maintains  angle  of  attack  and  elevator  deflection  to 

within  acceptable  limits  for  an  initial  velocity  perturbation  of  20  ms 

is  p=it/900  rad  m  *s.  Equation  (87)  was  again  used  to  calculate  r^  for  the 

values  of  p  shown  in  figure  26,  excepting  again  for  rt/450  rad  m  ^s,  and 

values  of  v  were  selected  to  be  0,  0.2,  0.4.  The  minimum  value  for  r 

occurred  for  V=() .  2  and  (i-U/900  rad  m  V,,  for  which  t  tie  corresponding 

weights  are  V,  =10  ms  ,  ill,  =10  degrees,  and  a,  -2  degrees.  Trajectories 

obtained  for  this  tighter  control  of  velocity  are  compared  with  the 

trajectories  obtained  for  light  attitude  control  in  figures  29,  30  and  11. 

Further  comparisons  of  the  response  of  these  two  kinds  of  control  are  shown 

in  figures  32,  3 3  and  34,  where  angle  of  attack,  velocity,  and  altitude  are 

plotted  for  the  first  20  s  of  flight,  following  an  initial  velocity 

_  1 

perturbation  of  20  ms  '.  Figures  32  and  33  show  now  velocity  is  decreased 
by  Increased  incidence  but,  from  figure  34,  at  the  expense  of  increased 
attitude  angle.  The  lower  damping  of  the  short  period  mode  for  the  tightly 
controlled  attitude  case  is  evident  from  figure  32. 

The  tightly  controlled  velocity  case  is  more  appropriate  for  the  gliding 
projectile  because  it  suppresses  attitude  errors  and  provides  some  altitude 
correction.  Only  three  stall’s  are  fed  back  to  the  elevators,  these  arc 
velocity,  attitude,  and  angle  of  attack.  The  values  of  derived  feedback 
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gains  are:** 


-7.99  x  10  4  radian  m  *s 

4.39 

-4.39 


From  the  computed  trajectories,  the  errors  at  25  km  ground  range  resulting 
from  an  initial  velocity  perturbation  of  20  m  s  *  are  270  m  in  altitude  and 
0.12  degrees  in  attitude.  These  errors  require  a  full  beam  angle  of  5.6 
degrees  for  seeing  the  target  or  0.28  degree  per  in  s  *  of  initial  velocity 
error.  Maximum  attitude  error  is  0.85  degree  or  0.042  degree  per  in  s  '  of 
initial  velocity  error.  Improvement  in  altitude  error  is  possible  only  at 
the  expense  of  further  undamping  of  the  short  period  mode  but  tnis  is 
undesirable.  Stability  was  further  checked  by  computing  trajectories  for 
flight  with  an  initial  velocity  error  of  20  ms1  and  with  a  series  of 
random  wind  gusts  having  magnitudes  in  the  range  0  to  10  in  s  These 
gusts  are  shown  in  figure  35  and  the  responses  of  angle  of  attack,  attitude 
and  altitude  are  shown  in  figures  36,  37  and  38.  Figure  36  shows  the 
relatively  low  damping  of  the  short  period  mode  as  revealed  by  the  angle  of 
attack  response  to  the  first  few  random  wind  gusts.  For  these  gusts,  the 
angle  of  attack  perturbation  exceeds  the  maximum  value  limit  of  2  degrees, 
reaching  a  maximum  of  5.7  degrees  at  13.5  s.  This  is  not  surprising  since 
the  system  was  designed  to  control  velocity  errors  up  to  only  20  m  s  1 
within  angle  of  attack  perturbations  of  2  degrees.  As  in  the  case  of  the 
full  state  feedback,  failure  of  attitude  feedback,  undamps  the  short  period 
mode,  but  is  otherwise  stable  in  the  face  of  velocity,  or  angle  of  attack 
feedback  fa  lure. 


This  reduced  order  approach  to  optimal  control  system  design  resulted  in 
relatively  simple  algebraic  relationships  between  closed  loop  nole 
positions,  optimal  feedback  gains,  and  state  and  control  variable  weights. 
The  simpler  model  does  provide  insight  into  the  physical  behaviour  of  the 
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el  eva  lor- control led  glider,  indicating  the  difficulties  in  obtaining 
altitude,  attitude,  and  velocity  regulation,  simultaneously.  Obviously  the 
approach  cannot  safely  be  used  in  total  isolation  from  the  complete  system 
as  the  discovery  of  short  period  mode  instability  has  shown. 

/  Single  State  feedback 

Because,  as  mentioned  previously,  there  may  be  difficulty  in  setting 

initial  attitude  reference  axes  accurately  enough,  making  in-flight 

attitude  measurements  too  inaccurate  for  regulating  gliding  flight,  it  is 

of  interest  to  discover  what  can  be  achieved  by  velocity  feedback  alone. 

In  this  section,  the  response  of  the  complete  system  to  initial  velocity 

errors  is  considered  when  only  velocity  is  fed  back  to  the  elevators.  The 

search  for  a  velocity  feedback  gain  is  first  commenced  through  use  of  a 

modified  form  of  equation  <AS)  to  compute  Riccati  elements.  Using  the  same 

state  variable  weights  as  used  in  the  final  full  state  feedback  system  II 

design,  a  measure  of  the  quadratic  integral  performance  index,  excluding 

the  control  contribution,  is  obtained  by  computing  r)(  from  equation  (OS) 

with  the  quadratic  term  omitted  and  with  a.,j  replaced  by  a.^-b^K^  and  a([ 

by  a  |  t'b  jK  ^ -  I»  figure  39.  r^  so  obtained  is  plotted  against  K  j  and  a 

minimum  is  seen  to  occur  at  K^--6xl0  radian  m  s.  When  K j  is  negative 

the  elevators  reduce  angle  of  attack  for  positive  velocity  perturbations. 

As  discovered  previously  this  type  of  response  leads  to  the  best  altitude 

correct  ion  but  can  lead  to  irrecoverable  downward  plunging  flight  if  K  ^  is 

not  carefully  selected.  This  tendency  is  reflected  in  figure  39  where  a 

_  ^  | 

rapid  rise  in  t ■  is  shown  for  approaching  -/xiO  "  radian  m  s. 

However,  for  K^--t>xlO  *  radian  m  *r>.  the  short  period  mode  damping  was  too 
low  and  the  computed  trajectory  showed  instability.  Backing  off  from  this 
minimum  value  of  K ^  to  a  value  of  -5x10  i  radian  m  Js  produced  a  stable 
trajectory  with  a  low-damped  oscillatory  response  in  the  phugoid  mode. 
Plots  of  velocity,  attitude,  altitude  and  angle  of  attack  response  to  an 
init  i<.l  velocity  perturbation  of  20  m  s  1  are  shown  in  figures  40,  41.  42 


!>4  - 


and  4:).  These  figures  show  a  low-damped  oscillatory  response  having  an 
amplitude  in  attitude  perturbation  of  about  2  degrees,  and  in  altitude 
perturbation  of  about  270  in.  Maximum  angle  of  attack  perturbation  is  about 
1.8  degrees.  Figures  40  and  43  show  how  angle  of  attack  decreases  for 
positive  velocity  perturbation.  It  is  of  interest  to  compare  these 
trajectory  results  with  trajectores  obtained  for  a  positive  value  of 
when  angle  of  attack  increases  with  positive  velocity  perturbation.  For 
the  purpose  of  comparison,  trajectories  for  Kj=!>xl0  ^  radian  m  's  are  also 
shown  in  figures  40  to  43.  These  trajectories  show  that  amplitudes  ol 
initial  phugoid  response  are  increased  but  so  is  the  damping.  The  settling 
time  for  this  sort  of  control  is  small  enough  to  effectively  remove  initial 
perturbations  in  velocity  and  attitude  by  the  end  of  midcourse  flight,  hut 
its  accentuated  attitude  response  to  in-flight  velocity  disturbances 
renders  it  unsuitable  as  a  mideour.se  controller. 

The  controller  with  K^-IixlO  *  radian  m  *s  has  attitude  and  altitude 
responses  to  velocity  perturbations  which  have  smaller  values  than  the 
uncontrolled  system  ami  hence  this  sort  of  controller  could  find 
application  in  mid-course  flight.  However.  the  feedback  gain  for  this 
controller  is  not  fui  Iv  optimal,  having  been  selected  to  minimise  the 
response  to  initial  velocity  perturbations  only.  Nevertheless,  a  result 
given  in  reference  10,  could  find  application  in  providing  information  on 
how  such  a  sub  optimal  system  would  respond  to  a  general  set  of  initial 
conditions.  The  result  referred  to,  uses  the  fact  that  the  ratio  of  two 
quadratic  forms,  each  having  the  same  number  of  variables,  has  turning 
points  for  finite  values  of  the  variables  provided  that  none  of  the  square 
terms  disappear  from  the  denominator.  Thus,  for  the  ratio  x'tlx/x'Qx,  where 
Q  must  he  positive  definite,  it  is  shown  in  reference  10  that  the  maximum 
value  of  this  ratio  is  the  maximum  eigenvalue  of  Q  ^R.  To  see  how  this 
result  can  be  applied  to  a  control  system  we  return  to  equations  (.39)  and 
(40)  and  exclude  wind  disturbances  and  the  final  state  weighting  matrix  S. 
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Put  ting  n  |  ■=  —  K  x,  whore  K  is  now  any  gain  vector  which  does  not  deslahi  I  i  :;< 
tin*  closed  loop  system,  we  have 


x  =  AjX 


and  Jj(x,t)  = 


'x  (l)Q1x(l)dl 


=  x  (t)RjX(t) , 


where  =  A-bK 


Q1  =  Q+KPK 


and  is  constant 


Following  reference  10, 


(dJJ/dt)/J1  =  -(xTQ1x)/(xTR]x) 


where  T  is  the  maximum  eigenvalue  of  *R  .  Integrating  both  sides  of 
inequality  (90)  leads  to  the  inequality 


t  S  -T  loge{J1(x,t)/J1(x0)} 


The  usefulness  of  this  result  depends  on  the  value  of  the  upper  bound  T 
which  has  been  derived.  It  should  be  noted  that  the  choice  of  for 
determining  T  is  arbitrary,  apart  from  the  requirement  that  be  positive 
definite.  When  Qj  is  positive  non-negative,  the  maximum  eigenvalue  of 
Q  becomes  infinite,  rendering  (91)  a  useless  piece  of  information. 
From  equations  (88)  and  (89),  R,  corresponding  to  arbitrary  Q  can  be  found 


from  l  lie  matrix  equation 


R^Rj+RjAj+Qj  =  0.  (92) 

Suppose  now  that  we  have  a  system  with  known  closed  loop  behaviour,  perhaps 

an  optimal  system  with  full  state  feedback.  Let  the  known  value  for  the 

performance  index  of  this  sytem  at  time  tf  be  J(x(tt)).  Now  change  the 

feedback  gain  vector  K  to  another  set  of  values  such  that  closed  loop 

stability  is  maintained.  So  doing,  establishes  .  Suppose  further,  thaL 

a  reasonable  value  of  T  has  been  found  by  selecting  arbitrary  and 

solving  equation  (92)  for  R  ,  then  inequality  (91)  can  be  used  as  a  basis 

for  finding  an  upper  bound  for  the  time  the  changed  system  would  take  to 

bring  the  performance  index  J(x(t))  to  some  fixed  value  J(x(tf).  To  use 

inequality  (91)  we  need  to  fci.  l  the  maximum  value  of  J^(x,t)  which  just 

fits  inside  J(x(tf)).  The  mathematical  problem  then  is  to  find  the  maximum 
T 

for  J^(x,(t))~x  RjX  subject  to  the  constraint 

•J(x,(t))  -  xTRx 

=  J(x(tf)) ,  a  *  i  ant. 

This  is  achieved  by  introduce •  ;nge  multiplier  A  and  seeking  the 

turning  points  of  the  auxiliary  ium,uon 

J,(x(t))  -  \{J (x(t) )  -  J(x(tf))i 
which  occur  when 

d.)  /3x->aj/ax=o, 

or,  that  is  when  2RjX-2ARx-0,  (93) 

which  shows  that  A  is  an  eigenvalue  of  R  . 

T 

fre-mul tip lying  equation  (93)  throughout  by  x  ,  shows  that 
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J 1 -AJ=0 , 

and  tlms  it  J  is  the  maximum  value  of  J  which  just  fits  inside  J(x(tg.)) 

.1  =AJ(x(t  )) 

1m  t 

where  A  is  the  minimum  eigenvalue  of  R  *Rj, 

and  from  inequality  (91), 

tH  loge{XJ(x(tf))/J1(x0)}  (94) 

where  t  is  the  time  taken  for  the  modified  system  to  bring  the  performance 
index  J  of  the  initial  system  from  J(x(0))  to  J(x(t^)).  Inequality  (94) 
thus  affords  some  means  of  comparison  between  settling  time  for  a  given 
optimal  control  system  with  complete  state  feedback  and  the  stable  system 
with  the  same  plant  matrix  but  with  arbitrarily  selected  feedback  gains. 

8.  DISCUSSION  AND  RESULTS 

Linear  optimal  control  theory  has  been  used  to  determine  the  state  variable 
feedback  which  best  regulates  mid  course  flight  of  a  gliding  aircraft-like 
projectile  with  elevator  controls.  A  set  of  differential  equations  in 
perturbations  of  the  state  and  control  variables  about  the  desired  operating 
conditions  for  the  projectile  have  been  derived.  These  equations  provided  the 
linear  mathematical  model  necessary  for  application  of  optimal  control  theory 
to  determining  feedback  gains.  The  desired  operating  condition  was  the 
trajectory  resulting  from  the  constant  elevator  setting  for  which  the 
horizontal,  range  over  an  altitude  drop  from  4  to  1  km  was  maximised. 

in  gliding  trajectories,  aerodynamic  lift  and  weight  are  approximately  equal 
and  opposite  and  characteristically,  conditions  change  slowly  along  the  flight, 
path.  Consequently  it.  was  discovered  that  a  time- invariant  system  model  was 


-  58  - 


adequate  for  feedback  control  design.  Because  lift  depends  on  the  square  of 
velocity,  perturbations  in  velocity  can  easily  destroy  the  balance  between 
lift  and  weight  causing  the  projectile  to  rise  or  fall  about  the  desired 
flight  path.  To  control  those  excursions  about  the  desired  flight  path,  the 
elevator  controls  modify  flic  aerodynamic  lift  by  causing  changes  in  angle  of 
attack.  Because  at  a  fixed  horizontal  range  from  its  point  of  launch,  the 
projectile's  seeker  beam  must  be  able  to  capture  the  target  it  is  also 
important  to  restrict  angular  excursions  in  attitude  of  the  beam  axis  about 
its  direction  on  the  desired  trajectory. 

The  physics  of  regulating  altitude  and  attitude  in  response  to  velocity 
perturbations  through  use  of  feedback  to  elevator  controls  has  been 
considered.  It  has  been  discovered  that  there  are  three  extreme  ways  in  which 
elevator  controls  can  react.  Firstly,  velocity  can  be  controlled  by  bringing 
about  changes  in  altitude.  For  example,  if  velocity  is  too  high,  the  elevator 
controls  could  cause  the  projectile  to  climb  until  velocity  falls  to  the  value 
where  lift  just  balances  the  weight  and  flight  would  then  proceed  at  an 
altitude  consistently  higher  than  the  desired  value,  but  desired  values  of 
attitude  and  velocity  would  be  maintained.  The  second  way  in  which  elevators 
call  react  is  to  cause  a  decrease  in  augle  of  attack  in  response  to  increased 
velocity  and  vice  versa.  This  action  reduces  the  increase  in  lift  due  to 
increase  in  velocity  and  hence  tends  to  reduce  displacements  normal  to  the 
desired  flight  path  and,  at  the  same  time,  also  tends  to  maintain  the  desired 
attitude.  Reducing  angle  of  attack  in  response  to  higher  velocity  and  vice 
versa  leads  to  a  lower  mean  drag  than  that  for  the  uncontrolled  flight  at 
constant  angle  of  attack.  The  resulting  motion  is  oscillatory  ar.d  the  lower 
mean  drag  means  that  damping  is  lower  than  in  uncontrolled  flight.  If  this 
type  of  control  is  overdone  so  that  the  lift,  in  response  to  increased 
velocity,  becomes  less  than  the  weight,  then  the  flight  path  will  take  a 
downward  plunge  attempting  to  perform  a  180  degrees  backward  turn.  This 
behaviour  could  not  be  predicted  by  the  linearised  mathematical  model.  The 
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third  typo  of  rtevato  control  occurs  when  the  angle  of  attack  is  caused  to 
increase  in  response  to  increased  velocity  and  vice  versa.  This  results  in 
increased  mean  drag  which  attempts  to  reduce  velocity  to  the  desired  value. 
The  increased  differential  between  lift,  and  weight  initially  tends  to  increase 
displacements  normal  to  the  desired  flight  path  and  the  resulting  motion  is 
highly  oscillatory  with  greatei  amplitude  and  damping  than  in  uncontrolled 
flight.  The  oscillatory  nature  of  the  response  can  preserve  some  accuracy  in 
a  1 1 i tude . 

These  observations  on  the  three  extreme  types  of  response  lead  to  the 
following  conclusions.  Of  the  three  states,  velocity,  attitude,  and  altitude, 
it  is  possible  to  regulate  any  two  at  the  expense  of  the  third  using  elevator 
controls,  but  not  all  three  at  once.  The  control  system  design  method  used, 
effectively  mixes  the  three  different  types  of  response  in  proportions 
dependent  upon  the  weights  selected  for  the  state  variable  perturbations  in 
the  integral  quadratic  performance  index.  It  is  not  possible  to  obtain 
continuously  improving  flight  regulation  from  feedback  gains  obtained  by 
simultaneously  increasing  the  weighting  of  all  state  perturbations  in  the 
performance  index.  Because  altitude  and  attitude  are  the  important  parameters 
in  meeting  accuracy  requirements  of  the  mid  course  trajectory,  velocity  was 
not  heavily  weighted. 

Three  control  systems  have  been  considered.  In  the  first  system  all  five 
states  are  fed  back  to  the  elevators.  in  the  second  system,  based  on  a 
reduced  order  mathematical  model,  only  three  states  are  fed  back.  Thirdly,  a 
system  which  fed  back  only  velocity  to  the  elevators  was  considered  and  in 
this  case  the  performance  index  was  based  on  velocity  perturbations  only. 
Features  of  the  three  sytem  are  given  in  Table  2  below. 

Comparative  performance  of  these  systems  is  given  i .:  Table  3  below 
In  these  three  systems,  elevator  deflection  an;  ; e.s  are  less  than  0.7  radian 
and  angle  of  attack  variations  do  not  exceed  2  cegrees  in  response  to  velocity 
perturbations  of  up  to  20  ms  *.  All  three  systems  were  found  to  be  stable  in 
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the  face  of  successive  random  wind  gusts  having  magnitudes  in  the  range  0  to 
10  m  s  * .  It  was  found  that  time-invariant  models  of  system  behaviour  were 
adequate  for  determination  of  feedback.  Consequently,  the  full  state  feedback 


TABLE  3.  COMPARATIVE  PERFORMANCE  OF  3  FEEDBACK  CONTROL  SYSTEMS 


System 

Maximum  amplitude  of 

response  to  unit 

velocity  perturbation 

Amplitude  of  response 

at  end  of  midcourse 

flight  to  unit 

initial  velocity 

perturbation 

Total  seeker 

beam  angle 

required  per 

unit  initial 

velocity 

Alt  itude 

At  t i tude 

Altitude 

Atti  tude 

perturbation 

1 )  Full  S La te 

feedback  (11) 

lOm/ms  1 

0.0S°/ms~1 

lOm/ms'1 

0.022°/nis‘1 

0 . 24° /ms  * 

2)  Reduced 

order  model 

13.5m/ms  * 

0.04° /ms-1 

13.5m/ ms  1 

0.006o/ms_1 

0 . 280/ms" 1 

3)  Velocity 

feedback  alone 

20m/ms  ' 

0. 10°/ms  * 

15m/  ms 

0.063°/ms  1 

0.43‘Vms'1 

4)  No  feedback 

24m/ms  * 

0.37°/ms_1 

15m/ ms  1 

0.086°/ms*1 

0.47°/ms'1 

_ 1 

system  could  be  synthesised  bv  compensating  networks,  designed  by  classical 
methods,  when  not  all  states  are  available  for  measurement.  The  system  based 
on  a  reduced  order  model  performs  very  well  in  suppressing  attitude  response 
to  velocity  perturbations  and,  from  Table  3,  this  is  reflected  in  the  seeker 
beam  angle  size  needed.  However,  Table  2  shows  the  short  period  mode  to  be 
very  lightly  damped  and  the  persistence  of  short  period  oscillations  may  be 
undesirable.  The  velocity  alone  feedback  system  achieves  a  moderate 
suppression  of  attitude  response  but  is  not  effective  in  altitude  control. 
This  is  because  the  altitude  control  being  sought  is  control  with  respect  to 
horizontal  range  and  not  time.  With  the  full  feedback  control  system,  the 
correction  of  altitude  as  a  function  of  time  is  excellent  for  the  trajectory 
with  an  initial  velocity  perturbation.  However,  as  previously  pointed  out, 
for  elevator  controls,  regulation  of  altitude  and  attitude  with  time  can  only 
be  carried  out  when  velocity  regulation  is  slackened.  Thus  the  rate  of 


descent  may  be  excellently  regulated  but  the  perturbation  in  the  horizontal 
component  of  the  initial  velocity  will  tend  to  persist,  leading  to  a 
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persistent  altitude  error  with  horizontal  range. 

The  full  state  feedback  system  induces  well  damped  responses  to  disturbances, 
and  effectively  halves  the  seeker  beam  angle  that  an  uncontrolled  projectile 
would  require.  The  performance  of  this  sytem  indicates  the  measure  of  error 
suppression  obtainable  when  angle  of  attack  variations  are  limited  to  2 
degrees  about  the  operating  value  and  elevator  deflection  is  limited  to  0.7 
radian.  When  initial  trajectory  errors  cannot  be  confined  within  acceptable 
limits  a  different  strategy  in  selecting  a  mean  flight  path  is  called  for.  To 
meet  a  wide  range  of  initial  trajectory  errors  it  may  be  necessary  to  arrange 
automatic  on-board  selection  of  compensating  fixed  elevator  settings.  The 
penalty  for  this,  is  that  in  addition  to  the  requirements  of  the  feedback 
control  system,  more  on-board  computer  memory  is  needed  together  with  elevator 
setting  algorithms. 
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LIST  OK  SYMHOI.S 


A 

It 


"DO 


m 


F 


o 


H 

i 


K 


K 


t 


L 


I. 

c 


M 

V 

Q 

i< 


s 


T 

s/(nr.AS  ) 

W 

V 

V 
X 


I’ I  ant  matrix 

Moment  of  inertia  of  projectile  in  pitching  plane  aliont  axis 
through  centre  ol  gravity 
Drag  coefficient 

Drag  coefficient  at  zero  angle  of  attack 

Lift  coefficient 

Ditching  moment  coefficient 

Aerodynamic  drag,  or  in  equation  a  disturbance  matrix 

Unmeasurable  state  selection  matrix 
Unconstrained  fund ional 

Scale  height  for  exponential  atmospheric  density  approximation 
Unit  matrix 

Integral  quadrat  ie  performance  index 
feedback  gain  vector 

Feedback  gain  vector  used  in  estimation  process  to  minimise  (x-x ) 

Lift,  force  on  projectile 

Lift  contribution  from  controls 

Pitching  moment  about  centre  of  gravity  of  projectile 
Control  variables  weighting  matrix 
State  variables  weighting  matrix 
Kiccati  matrix 

Reference  area  for  aerodynamic  coefficients,  or  in  equation  (40), 
a  final  state  weighting  matrix 

Upper  hound  ol  settling  time  for  a  dynamical  system 

Induced  drag  factor 

Velocity  vector 

Velocity  magnitude 

Horizontal  range  variable 


a  . 


F.lement  ol  the  plant  matrix  A 


Control  vector 


Component k  of  cont  rol  vector  h 

Components  of  control  vector  for  reduced  order  system 

-  ■> 

Gravitiition.il  constant,  9.8  in  s 
Alt  it inlc 

Unit  vector  in  horizontal  downrangc  direction 
Unit  vector  in  horizontal  direction  normal  to  i 

e 

Unit  vector  direction  vertically  downwards 
Velocity  feedback  gain  in  reduced  order  system 
Flight  path  angle  feedback  gain  in  reduced  order  system 
Integer 

Mass  of  projectile 

Unit  vector  normal  to  flight,  path  forming  right  Hand  set  with  f 
Angular  pitching  rate 
K lenient  of  R i ccat  i  matrix  R 
T  i  me 

Unit  vector  in  a  forwards  direction  tangential  to  flight  path 

Measurement  noise  vector 

System  noise  vector 

State  variable  vector 

Ks t  imate  of  x 

State  variable  vector  for  reduced  order  system 

Angle  between  longitudinal  axis  of  glider  and  horizontal 

Angle  of  attack 

Elevator  deflection  angle 

Perturbat  ion  in  (•) 

Variable  used  in  characteristic  polynomials  or  Lagrange  multipl 
Relative  velocity  weight  for  reduced  order  system 
Angle  between  atmospheric  wind  and  the  horizontal 
Relative  flight  path  angle  weight  tor  reduced  order  system 


Atmospheric  density 


(»  I  I'nnslitiil  reference  atmospheric  density  corresponding  to 

Stale  lt«' i  kIK  II 

o  Him  I  part  of  complex  root  of  characteristic  cipiation 

to  Imaginary  part  of  complex  root  of  characteristic  equation 

T  Dummy  t  fine  va r  i  ah  l e 

»|i  flight,  path  direction  angle  ((-)*<») 


Superscripts 


T  Transposed  matrix 

Suhscr i pt s 

0  Unperturbed  01  optimal  trajectory  conditions 

I  Perturbed  trajectory  conditions 

i,j,k  Integers 

in  Denotes  maximum  value  of  a  variable  for  weight  selection  purposes 

■  Initial  va i uo  a  i  t  =0 

f  At  I  i  na  I  va  I  tie  of  t  ime 

w  Wind 

liw  Horizontal  wind 

vw  Vertical  wind 

K  He  la live  wind 

V  Velocity 

oi  Angle  ol‘  attack 

0  Attitude 

t|  Pitching  rate 


h 


Alt i tnde 
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APPENDIX  I 


PROJECTILE  MASS  AND  AERODYNAMIC  PROPERTIES 


Mass 


=  280  kg 


Pitching  moment  of  inertia  B 


=  73  kg  m2 


Reference  area 


=  0.06  m2 


Reference  length 


=  0.28  m 


Wing  aspect  ratio 


=  7 


Wing  area 


=  0.9  m' 


Induced  drag  factor  S/jtkAS  =  1.378  x  10 


9.8ms 


-1 


Sea  level  air  density  p 

ref 


=  1.228  kg  m 


Air  density  at  altitude  h,  p 


Pre£  exp(-h/M),  where  H=ll)4m 


p  ,S/2m  -  1.3157  x  10"4  m"1 

r  ref 


PrefSb/2B  -  1 -A130  x  10'4  m'2 


Pr(>fSb2/4B  =  1.9783  x  I0_r>  nf1 


Aerodynamic  Force  Coefficients 


(p  S/2m)C  =  3.2893  x  10"5  m'1 

aCj/3a  =  25 

(p  S/2m)3Cr/3a  =  3.2893  x  10‘3  nT1 

3Cr/3n  -  1 

(S/neASw)(3C[/3a)2  =  8.612 

C  (a)  =  Cno+S/neASw  ( 3CL/3a)2  a2 
3CD.  =  (2S/ncASw)(3Cj /3a)2  a 

Aerodynamic  Moment  Coefficients 

3C  /3a  =  -13 

m 

(p  ,Sb/2B)3C  /3a  =  -1.837  x  10‘J  m*2 

rref  m 

3C.  /3q  =  -112 

m  1 

(p  -Sb2/4B)3C  /3q  -  -2.2156  x  H)"3  m” ' 

'  rrf  m 

3C  /3a 

it< 
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(prefS,,3f/4B)acIn/ao  =  '4  x  lo'4  m" 1 

8cB/3n  =  -4 

(prefSb/2B)  ^V9'1  =  -5- 652  x  lo'4  m"2 

Optimal  Glide  Path  Par  ;rs 

From  equation  (31),  ^  =  0.17037  rad  =  9.76  degrees 

and  from  equation  (27),  Vfl  =  141.14  exp  (h/2H)  m  a'1, 

with  cos (0  -  a  )~1.  rT  , 

Maximum  value  of  Lilt/Drag  =  »  (ac  /{*») 

o  L 

C0o*“o(8cJ  ./aO^/neAS, 

♦  n0OcL/on) 

CDo+°^8Ci/8°)2s/«ASw 

=  7.411  tor  «o  =  0.  17037  rad 
and  r)0  =  -0.5537  rad. 

From  equation  (29),  tan^-of^)  =  - 1 . 3224/ { (V2/2H)  }  (j  2) 

From  equation  (27),  range  X  --  7.411  {(V  2-V  2/19  6) 

°  oi  o  '  '  ' 

+(hoi-ho)i  m  (1.3) 

and  from  equation  (30),  time  t 

=  7.41I[{(Vo.-Vo)/9.8}+2H{(1/V()-(1/Vo.)}] 


(1.4) 


Also,  sujuaiion  (2R)  gives  V 


=  (V  /2H)s  i ii (0  "Of  ) 
o  o  o 


(1.5) 


i\  =  4  x  10J  m 

i 


103  m 


V  .  --  I  72.  4  m  s 

01 


V  t  :  148.4  m  s 

ot 


0  -of  .  -  -0.11664  rad 

01  oi 


0  ,-Ct  _  =  -0.12072  rad 

of  of 


=  0.05373  rad 


0  ,  =  0.04965  rad 

of 


XQf  =  25,144  m 


tf  =  157.2  s 


=  -0.1729  m  s 


V  ,  =  -C.1326  m  s 

of 


Values  of  elements  a_  of  plant  matrix  A  as  given  by  equations  (34)  to  (39) 
for  initial  flight  conditions  (h^=4xl03m) 


1  -0.0152s  2.041ms 


-9.733ms 


0.000131  Is 


2  -0.00065429m"1  -0.38673s"1  0.0066153s"1  1 


0. 000006m" 's-1 


3  0 


0 


0 


1 


0 


■"'#8 
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APPENDIX  I t 

CIIAKACTKKISTIC  EQUATION  AND  SOLUTION 

In  this  append i x  the  characteristic:  equation  corresponding  to  equation  (39), 
with  coefficients  assumed  to  he  time-invariant,  is  solved  by  an  approximate 
method  so  that  the  dynamical  behaviour  of  l.he  glider  can  be  related  to  its 
aerodynamic  and  mass  properties.  By  comparing  the  magnitude  of  the  most 
significant  elements  a.,  of  the  plant  matrix  A,  given  in  Appendix  l  for 
initial  and  final  values  of  flight,  it  can  be  seen  that  the  assumption  of  time 
invariance  is  reasonable.  This  situation  is  brought  about  because  the 
aerodynamic  coefficients  have  been  assumed  to  have  constant  values  over  the 
speed  range  and  in  gliding  (light  the  dynamic  pressure,  given  by  pVa/2,  is 
then  approximately  constant.  The  Mach  number  range  based  on  optimal  glide 
path  conditions  varies  from  0.53  initially  to  0.44  finally,  which  is  well 
outside  the  transonic  range,  so  that  the  assumption  of  constant  aerodynamic 
coefficients  is  reasonble. 

It  is  evident  from  the  very  small  values  of  the  elements  in  the  fifth  column 
of  the  plant  matrices  given  in  Appendix  I  that  the  state  variables  are  only 
weakly  dependent  upon  altitude.  Hence  for  open  loop  analysis,  at  least,  it  is 
reasonable  to  consider  only  the  fourth  order  system  obtained  by  taking  the 
elements  in  the  fifth  column  of  A  to  be  zero.  If  the  roots  of  this  fourth 
order  equation  are  A^,  ,  A.^,  A^,  then  the  fifth  root  corresponding  to  the 

complete  quint ic  is  approximated  by 

\  =  IAI/AjAA^  (11.1) 

where  |A|  is  the  determinant  of  the  complete  A  matrix.  In  considering  the 
system  given  by  equation  (39)  it  is  worth  including  feedback  of  all  the  state 
variables  excluding  altitude.  Thus  equation  (39)  will  be  treated  as  the 


system 


X  ~ 


|k1k2k3k4  |x 


(11.2) 
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al  1 

a  1 2 

a  1 3 

a  14 

r  bl 
1 

• 

• 

x  - 

b2 

• 

b3 

1 - 

z> 

a42 

a43 

a44 

_ 

b4 

or  x  -  Cx 


(11.3) 


whore  the  elements  of  the  matrix  C  ate  given  by 

c.  .  -  a.  .  -  l> .  K  . .  (11.4) 

ij  t.)  i  .1 

The  characteristic  polynomial  corresponding  to  the  system  given  by  equation 
(11.3)  is  found  by  expanding  the  determinant 

fo  no  | A I 4  -  C|  in  the  form 


X-  ♦  C„A»  *  C,A*  ♦  C2A  »  C3. 


where  it  is  found  that 


C0  =  ‘(all+a22+a44)  +b2K2+b4K4 


(11.5) 


C1  =  all(a22+a44)  ’a22a44'a24a42_al2a2l'a43 


+b2Kla12'K2(allb2+a44b2'a24b4)  +b4K13 


■K4(allVa22b4'a42b2) 


(11.6) 


C2  ~  a11a24a42'alla22a42+al2a21a44"aI2a41a24 


*a42a22"a42a23 
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~K1(a12a44b2~a12a24b4~a13V 

”S(alla24b4“a11a44b2+a43b2“l‘23b4) 

“K3Ullb4'a42b2+a22b4)  (lI*7) 

”K4(aila42b2~a1la22b4+a12221b4'ai2a41b2) 


C3  =  alla23a42*alla22a43+a12a21a43"al2a41a23 


+a13a4la22~a13a42a21 


(II. 8) 


'Kl(a12a43b2'a12a23Va13a22Va13a42b2) 


■K2(aHa23b4'aHa43b2+a13a41b2'a13a2.b4) 


‘K3(alla42b2“aUa22b4+a12a21b4'a12a4lb2) 


By  carrying  out  an  order  of  magnitude  analysis  of  these  expressions  for  the 
coefficients  from  the  values  of  a .  ^  given  in  Appendix  I  it  is  possible  to 
simplify  the  results  as  follows:- 


C0  =  "^a22+a44^+b2K2+b4K4 


(II. 9) 


C1  =  'a42+b2Kla12+b4K2+b4K3 


'K4*b4(aH+a22)"b2a42* 


(11.10) 


C2  =  a42(aira23)+Klb4(al3+a12)'K2b4(aira23) 


'K3'b4(all+a22)~a42b2* 


'K4*b4*a12a2ralla22^+b2alla42* 


(11.11) 


+K2b4a13a21 


'K3{b4(a12a2ralla22)+b2aHa42^ 
II. I  SOLUTION  OF  QUART I C 

An  approximate  solution  to  the  quartic 


(11.12) 


A4  *  C0A;l  t  C jA2  +  C2A  +  C3  =  0  (11.13) 

can  he  obtained  by  the  method  of  reference  13.  Using  the  notation  of 
reference  13,  the  quartic  can  be  written  as  the  two  quadratic  factors, 

(A2+(a-a)A+p-b)(A2t(a+a)A+p+b)=0  (11.14) 

where  c t  =  C^/2  >  P  =  (A+£)/2  , 

a  =  il/2  ,  b  =  a(«-B/£)/2  ,  (11.15) 

A  -  Cj  -  a2  ,  B  =  C2  -  «A 
ami  f,  satisfies  the  cubic  equation 

£3  +  (2A-a2)£2  +  (A2+2Ba-4C  )£  -  B2  =  0  (11*16) 


The  approximate  solution  of  this  equation  is  described  in  reference  (13). 
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By  substituting  £  =  x(-B)2^3  in  equation  (11.16)  it  is  found  that  the 

resulting  cubic  equation 


x3  +  (2A-«2)(-B)*2/3x2+(A2+2Bu-4C3)(-B)‘4/3x-I=0  (11.17) 

has  one  real  root  for  0<x<l.  Over  this  interval,  the  Chebyshev 
approximation  x'^-1 . 5x2-0 . 5625x+0 . 03125  is  substituted,  reducing  equation 
(11.17)  to  the  quadratic 

{(2A-a2)(-B)‘2/3+1.5}x2  +{ (A2+2Ba-4C3) (-B)"4/3-0 .5625 }x 
-1  +  0.03125  =  0,  (11.18) 


which  leads  to  the  solution  for  £, 


t  =  -  (A2+2Bflf-4C3)(-B)‘2/3-0.5625(-B)2/3 


2(2A-ttz)(-B)“2/3+3 

(A2+2Ba-4C3)(-B)'2/3-0.562.5(-B)2/3 

2(2A-«2)(-B)~2/3+3 


+  2(1-0. 03125) (-B) 


4/3 


2(2A-az)(-B)'2/3+3 


1/2 


(11.19) 


This  solution  can  be  greatly  simplified  when  the  following  conditions  hold, 


!(A2+2Ba«4C3)(-B)"4/3|  >>  0.5626. 

|(2A-«!l)(-B)“2/i|  »  1.5, 

|4B2(2A-«2)/(A2+2B«-4C3)2|  «  1,  (11.20) 

|2B«-4C3I  «  A2  . 

and  ignoring  0.03125  compared  with  1,  gives  the  simple  approximation 

i  =  B2/A2  or  iX/2  -  -B/A.  (11.21) 

In  summary,  an  approximate  set  of  solutions  to  the  qua rtic  equation  (11.13) 
is  given  by 

A.  .,  =  -(a-a)/2  ±  [(«-a)2/4  -  (p-b)]1/2, 

X3  u  =  -(«+a)/2  t  [ (a+a)2/4  -  (p+b)|l/2  (11.22) 

where  a,  p,  a,  b,  are  given  by  equations  (11.15)  and  the  parameter  £ 
approximated  by  equation  (11.19).  An  exceptional  case^  occurs  when  B--0 
in  equations  (11.15),  making  £=0  a  solution  of  equation  (11.16).  In  this 
case  the  limiting  value  of  b  is  0.5(A2-4C3)^2  unless  (A2-4C3)  is  negative 
in  which  case  the  positive  root  of  the  quadratic 

l.2  +  (2A“«2H  +  (A2-4C3)  -  0  (11.23) 

is  taken  as  the  value  of  £  in  equations  (11.15).  A  modification  to 

equations  (11.22)  should  he  used  when  p  and  b  are  approximately  of  the  same 
value.  In  this  case,  the  relationship 

p-b  =  C3/(p+b)  (IT .24) 

obtained  by  comparing  equations  (11.13)  and  (11.14),  should  be  used. 
Similarly  when  0  and  a  are  approximately  of  the  same  value,  the 


relationship 


a-,i  =  C^“20/(fl+a)  (11.25) 

should  be  used. 

Inserting  the  values  of  a..  from  Appendix  1  into  equations  (11.9)  to 
(11.12)  with  /.ero  feedback  gains,  it  is  found  that 

A-C, 

B  ~  C2  •  ClV2 

and  that  all  conditions  (11.20)  are  well  satisfied.  Hence,  from  equations 
(11.15)  and  (11.21) 

a  *  f,1/2  «  C0/2  -  C2/Cj, 
o+a  =  CQ  -  C2/Cj  ,  o-a  =  C^/Cj , 
b  =:  (l/2)jCl  +  (ty4)(C0-2<yC])}, 

P  =  (1/2)|C1+(1/A)(C0-2C2/CJ)!E} 

It.  is  found  that  ( Cy/2  -  tyCj)2  and  Cq/2  (Cg/2  -  tyCj)  are  negligible 
compared  with  ,  so  that 

P+b  “  Cj,  and  from  oqution  (11.24) 


p-b  =  C3/Cj 

From  equations  (II. 5)  to  (II. 12),  the  following  approximations,  in  terms  of 
the  elements  of  the  plant  matrix  A,  are  obtained 

a+a  =  •(a22+\4^a23)'  a~a  =  a23_a  1 1  ’ 


A I 


(HI,  t  -»42 


,  (Hi,  =  «,  V^r^.r'n 


(11.20) 


Substituting  these  results  into  i'i)u<ition  (11.22)  shows  the  short  period 
damping  and  frequency  to  be  approximated  by 

°sp  ~  (  l/2)(a22,i*4A+a2:i)  ’  ‘“"si,  =  "a42'<,2si>,(ll'2;) 


and  the  long  period  damping  and  frequency  hv 


C,lp  =  (l/2)Catl-a23),  =  5,i3a2ra23air0*lp- 


(11.28) 


Using  the  numerical  values  of  ;k  .  from  Appendix  I,  the  initial  open  loop 
values  of  damping,  frequency,  and  period  are, 

a  =  -:).455xl0'1s'1  ,  to  =  6.036  s'\  T  =  1.041  s, 
sp  sp  sp  * 


0(|)  =  -1  ,09lxlo"2s“1  ,  iu1()  =  7.969X10-V1  ,  Tj  =  78.85  s, 


and  the  final  values  are 


u  -  -4.015xl0"1s'1 ,  w  -  6.034  s'1,  T  =  i  .041  s, 
sp  ’  sp  sp 


O.  -  -1  .iAlxlo'^s  UI  =  8.634x10 '  's' 1  ,  T.  -  72.78  s. 

Ip  ’  sp  Ip 

The  elements  a.j  can  be  expressed  in  terms  of  the  significant,  aerodynamic 
coefficients  when  the  gravity  components  are  eliminated  through  equations 
(22)  and  (23).  This  results  in  the  following  approximations 


•<Vfl/2)  [  (pS/2,«)  { (8C,  /On)  +CD(«  )} 


♦  CpSh2/ 4B)  {  (0Cm/3qti)Cm/ 0« ) }  ] 


(11.29) 


iu  •  2:  - (pSb/2B)  V,,  (i)t:  /Dor) 

sp  0  m 


(II .30) 


-(3/2)  V0(pS/2m)  CD(aQ) 


(11.31) 


U»1J  ~  v'2  V^(pS./2m)  CL(a())  (II.  32) 

1 1 . 2  REDUCED  ORDER  SYSTEM 

The  prospect  of  a  reduced  order  system  to  represent  the  long  period  motion 
is  apparent  from  consideration  of  equations  (11.14)  and  (11.26).  The  first 
three  plant  equations  of  the  fourth  order  system  are 

Xj  =  a , | xj  +  a 12x2  +  a13x3 


a21Xl  +  a22X2  +  a23X3  +  X4 


and  these  can  be  written  as 
Z1  =  al lZl+a13Z?+(a12+a13)x2 


Z2  "  ■|,21zra23Z2"^a22+'*23^X2 


(11.33) 


where  Zj  =  x^ 

and  z ^  =  x^-  (11.34) 
Regarding  z^,  z^  as  controlled  variables  and  x0  as  the  control  variable, 
then  the  characteristic  equation  for  the  open  loop  part  of  this  system  is 


*  A(a23‘aU)  +  al3a2ra23ail  =  °’ 
which  from  equation  (11.26)  is  equivalent  to 


X2  +  A(oi-a)  +  p-b  -  0, 


the  long  period  quadratic  component  of  the  complete  quartic  given  by 
equation  (IT.  14).  In  terms  of  the  physical  variables, 


7. 


V 


the  velocity  perturbation, 


